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1.  RESEARCH  OVERVIEW 


We  have  continued  work  at  various  levels  of  our  IU  system  and 
started  to  evaluate  techniques  for  applications  to  DMA  supplied 
images.  This  report  constitutes  the  final  report  on  contract 
F-33615-76-C-1203.  The  last  six  months  work  can  be  described  under 
the  following  headings: 

Image  Matching 

Matching  of  an  image  to  a  symbolic  map,  or  a  symbolic  description 
of  another  image,  is  central  for  the  tasks  of  map  updating  and  change 
detection.  In  the  past  we  have  described  results  using  aerial  images 
of  areas  such  as  San  Francisco,  Stockton,  San  Diego,  etc.  fl].  These 
previous  techniques  used  a  simple  matching  scheme,  with  each  element 
in  one  description  being  matched  to  the  best  corresponding  element  in 
the  other  description,  and  not  allowing  for  any  revision  based  on  the 
matching  of  other  neighboring  elements.  We  have  now  incorporated  a 
relaxation  matching  algorithm.  This  algorithm  is  different  from  those 

used  by  Rosenfeld  and  associates  f2],  in  the  use  of  a  well  defined 
optimization  criterion. 

Texture  Analysis 

We  have  continued  development  of  our  structural  texture  analysis 
techniques.  The  analysis  uses  micro-edges  detected  in  a  texture  and 
derives  repetition  pattern  characteristics  of  these  edges.  Our 
previous  presentations  have  described  techniques  for  determining  the 
width  of  texture  primitives  and  a  repetition  period,  if  any.  Our  new 
techniques  are  also  able  to  extract  the  length  of  the  primitives  and 
thus  describe  the  shapes  of  the  primitives.  The  usefulness  of  these 
descriptions  for  recognition  of  natural  textures  is  currently  being 
tested . 
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Texture  Synthesis 


We  have  several  ongoing  projects  in  synthesis  of  natural  textures 
from  a  stochastic  model  using  few  parameters.  The  techniques  include 
auto-regressive  modeling  with  conditional  expectations  and  algebraic 
reconstruction  techniques.  Models  for  texture  synthesis  are  also 
useful  for  texture  analysis,  as  they  validate  the  sufficiency  of  the 
models  used  for  analysis. 

Segmentation 


We  are  trying  to  use  the  texture  analysis  techniques  to  aid  in 
scene  segmentation.  Texture  features  can  be  used  as  intensity 
features  for  segmentation.  However,  difficulties  arise  because 
texture  features  must  be  measured  over  a  window  assumed  to  contain  a 
single  texture  only.  Also,  texture  features  have  many  components  and 
should  be  treated  as  a  vector.  Faugeras  and  Lee  give  some  preliminary 
results  in  Section  2.4. 

Another  segmentation  project  is  to  develop  techniques  for 
segmenting  images  that  have  unimodal  intensity  (or  color)  histograms. 
This  happens  typically  when  the  image  consists  mostly  of  a  large 
background  region,  with  small  but  significant  other  regions.  Faugeras 
and  Bhanu  have  developed  a  gradient  relaxation  technique  to  modify  the 
histogram  to  bring  out  the  peaks  corresponding  to  the  small  regions 
(Section  2.9).  Currently,  this  technique  is  applicable  if  only  two 
types  of  regions  are  present  in  the  image. 

Range  Data  Processing 

Our  previous  linear  feature  extraction  techniques,  developed  for 
intensity  image  processing  (3),  have  also  proved  useful  for  range  data 
processing.  However,  for  range  data  we  are  able  to  obtain  complete  or 
near  complete  object  boundaries  by  subsequent  processing  that  exploits 
the  richer  information  contained  in  such  data.  The  DARPA  contract  has 


only  supplied  minimal  computing  funds  for  this  work.  However,  it  is 
included  here  because  of  its  potential  application  to  aerial  image 
processing . 


Hardware  Implementation 

In  continuing  work  with  Hughes  Research  Laboratories,  Malibu, 
California,  we  are  investigating  the  use  of  VLSI  technology  for 
hardware  implementation  of  IU  algorithms.  We  have  chosen  to 
investigate  the  following  algorithms  initially: 

i)  Nevatia-Babu  Line  Finder  [3] 

ii)  Ohlander  Region  Segmentor  [4] 

iii)  Laws  Texture  Analysis  System  f5] 

The  choice  of  the  above  three  algorithms  was  based  on  their 
computation  intensive  nature,  their  use  for  a  broad  range  of  problems 
and  experience  with  a  large  number  of  images  for  the  first  two.  Also 
these  algorithms  are  largely  local  and  hence  easier  to  implement  in 
VLSI  hardware,  where  reducing  interconnections  is  important.  Further, 
the  three  algorithms  have  common  kernels,  such  as  convolution,  but 
also  require  different  subsequent  processing.  A  study  of  there  should 
provide  valuable  feedback  on  the  feasibility  of  hardware 
implementation  for  a  large  class  of  algorithms. 

At  this  time,  no  decision  on  algorithms  for  actual  implementation 
has  been  made  and  opinions  of  the  IU  community  are  invited  on  the 
suitability  of  the  proposed  algorithms  as  well  as  suggestions  for 
other  algorithms. 
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2.  IMAGE  UNDERSTANDING  PROJECTS 


2.1  Semantic  Description  of  Aerial  Images  Using  Stochastic  Labeling 
O.D.  Faugeras  and  K.E.  Price 


Introduction 


The  purpose  of  computer  scene  analysis  is  to  automatically 
produce  a  description  of  the  content  of  an  image  similar  to  one 
obtained  from  a  skilled  human  observer.  In  order  to  achieve  such  a 
goal,  a  symbolic  description  of  the  raw  image  data  must  be 
constructed.  This  requires  the  application  of  many  of  the  now  well 
developed  techniques  of  Image  Processing  (image  bandwidth  compression, 
image  restoration  and  image  enhancement) ,  extraction  of  features  such 
as  texture  and  edges,  segmentation  of  the  image  into  homogeneous 
regions  with  respect  to  one  or  several  properties,  and  measuring 
features  that  characterize  these  segments  (color,  brightness,  texture, 
size,  and  shape)  and  also  relations  between  these  segments  (brighter 
than,  larger  than,  above,  below,  etc.).  The  output  of  this  complex 
sequence  of  processes  is  something  that  does  not  resemble  the  original 
input  array  of  pixels  but  is  much  more  suitable  for  high  level 
processing:  a  symbolic  description  which  is  represented  as  a  labelled 
graph  or  semantic  network.  To  proceed  any  further,  we  must  also 
assume  that  we  have  access  to  another  body  of  knowledge  containing  a 
priori  information  about  the  expected  content  of  images  of  a  given 
area.  We  will  not  make  any  assumptions  about  how  this  world  model  has 
been  obtained  (manual  input  or  intelligent  learning)  and  will  only 
assume  that  its  representation  is  the  same  as  the  image,  i.e.,  that  it 
is  also  a  semantic  network.  The  process  of  obtaining  a  semantic 
description  of  a  given  image  can  then  be  viewed  as  finding  the 
solution  of  a  graph  matching  problem:  either  match  the  image  onto  the 
model  or  match  the  model  onto  the  image.  Thus,  this  is  equivalent  to 
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the  general  graph/subgraph  isomorphism  problem  which  is  well  known  to 
be  NP-complete.  Practical  and  useful  solutions  can  nonetheless  be 
found  as  we  will  show  in  this  paper. 

Matching  techniques  must  be  described  in  terms  of  performance  on 
a  well  defined  problem.  The  particular  task  under  study  is  the 
analysis  of  an  image  of  a  scene  using  an  approximate  specification  of 
the  scene  which  would  apply  for  many  different  images  of  the  scene 
(i.e.  a  model).  Both  the  image  and  scene  are  represented  by  semantic 
networks,  the  image  description  is  automatically  generated  and  thus 
reflects  any  errors  in  the  segmentation  process.  The  model  is 
specified  by  the  user  and  contains  only  the  important  objects  and 
relations . 

A  similar  problem  was  attacked  by  Rubin  [1]  with  a  search 
procedure  and  a  more  detailed  model.  His  work  has  been  combined  with 
a  relaxation  procedure  by  Smith  [2],  The  general  form  of  the  solution 
is  similar  to  ours  but  the  scene  model  is  significantly  more  detailed 
so  that  exact  comparisons  cannot  be  made. 

We  will  first  present  the  image  and  model  descriptions,  which  are 
the  input  to  the  matching  procedure,  then  the  basic  matching  technique 
which  we  use  to  compare  two  objects.  Next,  the  global  matching 
process  (labeling)  will  be  described,  finally  we  will  discuss  the 
results  of  applying  this  procedure. 

Description  and  Basic  Matching  Techniques 

This  section  will  first  present  a  description  of  the  symbolic 
representation  used  for  the  scene  and  image.  Then  the  method  by 
computing  initial  likelihoods  of  particular  matches  is  described  along 
with  the  method  to  compute  compatibilities  of  pairs  of  matches. 
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Image  And  Model  Descriptions 


Our  matching  system  uses  a  feature  based,  symbolic,  description 
of  an  idealized  scene  (the  model)  and  of  the  image  of  a  portion  of 
this  scene  [3],  The  image  description  is  derived  automatical ly  from 
an  image  and  the  model  is  developed  by  the  user  through  an  interactive 
procedure . 

The  basic  objects  used  for  the  image  description  are  the  segments 
of  the  image  generated  both  by  a  general  region  based  image 
segmentation  procedure  [4]  and  by  a  linear  feature  extraction 
procedure  [5].  The  regions  are  derived  by  locating  connected  areas 
which  are  uniform  with  respect  to  some  feature  in  the  input  image 
(color  parameters,  texture  values  etc.)*  Linear  features  are  defined 
as  long  narrow  objects  which  differ  from  the  background  on  both  sides, 
and  are  described  by  as  a  sequence  of  straight  line  segments  with  some 
small  width.  Typically,  the  images  which  we  use  have  a  total  of 
100-200  individual  segments  of  both  types.  The  symbolic  description 
is  completed  by  computing  various  features  of  the  segments  and 
relations  between  them. 

The  features  used  for  the  symbolic  description  are  those  which 
can  be  reliably  computed  from  the  available  data  (the  input  image  and 
the  region  or  line  descriptions).  They  include  properties  such  as 
average  color  and  texture  (currently  only  simple  texture  measures), 
size,  position,  two-dimensional  orientation,  and  simple  shape 
measures.  Also  included  are  various  relations  between  image  segments 
such  as  adjacency,  nearby,  and  relative  positions  (above,  below, 
etc.).  With  all  these  relations  a  segment  may  easily  be  related  to  as 
many  as  100  other  segments.  This  description  is  not  intended  to  be 
used  for  reconstruction  of  the  original  image,  it  is  meant  to  capture 
the  important,  observable,  information  contained  in  the  image. 

The  model  description  is  identical  to  that  used  for  the 
image  -  feature  based  descriptions  of  basic  region-like  and  line-like 
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elements  including  relations  between  them.  Additionally,  the  basic 
elements  in  the  model  are  grouped  into  more  complex  objects, 
associated  with  generic  descriptions,  and  referred  to  by  actual  names. 
The  feature  values  in  the  model  will  not  correspond  exactly  to  image 
values,  but  are  approximations  of  the  likely  values,  so  that  one  model 
of  a  scene  can  be  used  with  many  similar  images  of  the  same  scene. 
The  relations  between  elements  which  are  included  in  the  model  are  the 
"important"  ones,  that  is,  if  a  relation  appears  in  the  model 
description  then  it  is  expected  to  occur  in  the  image  description, 
but,  if  no  relation  occurs  in  the  model  then  nothing  may  be  said  about 
its  appearance  in  the  image  (negative  relations  could  be  used,  such  as 
must  not  be  adjacent,  etc.).  Similarly,  only  the  important  objects 
are  described  in  the  model,  thus  it  is  not  a  complete  description  of 
the  entire  scene.  Generally,  the  model  description  is  smaller  than 
the  image  description  containing  20-30  basic  elements. 

In  summary,  the  input  for  the  matching  procedures  is  the  symbolic 
descriptions  of  both  the  input  image  and  a  mode)  of  the  scene.  The 
model  description  determines  the  outcome  of  the  matching  operations. 
The  image  description  is  automatically  derived  and  may  contain 
errors  -  especially  where  simple  objects  are  broken  into  several 
pieces.  The  model  description  is  incomplete,  as  it  contains  only 
important  objects,  thus  most  segments  in  the  image  (and  objects  in  the 
actual  scene)  will  not  be  described  by  objects  in  the  model. 

Basic  Matching  Technique 

The  global  matching  procedure  is  a  stochastic  labeling  procedure 
(also  called  relaxation  procedure)  which  will  be  explained  in  detail 
in  the  next  section.  The  relaxation  technique  requires  a  basic 
procedure  to  compare  how  well  a  model  element  agrees  with  an  image 
segment  for  its  operation.  In  our  previous  work  in  matching  pairs  of 
images  [2]  and  in  analysis  of  images  using  models  of  the  scene  Tfil  we 
have  developed  a  comparison  procedure  which  can  rate  the 
correspondence  of  an  object  in  one  image  (or  model)  with  an  object  in 
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another  image  (or  model).  The  basic  procedure  combines  differences  in 
all  available  features  and  relations  to  produce  a  single  rating  of  the 
quality  of  the  match.  The  past  experiments  indicated  that  this 
procedure  produces  reliable,  and  generally  accurate  measures  for  the 
differences  between  two  objects  (i.e.  the  mode]  based  matching 
performed  accurately  on  a  variety  of  scenes) .  The  problems  with  the 
past  matching  system  arose  in  the  use  of  these  results  by  the  global 
matching  procedure,  particularly  in  requirements  for  ordering  the 
selection  of  elements  to  match  and  the  handling  of  objects  which  break 
into  several  pieces.  (This  is  discussed  in  more  detail  later). 

Briefly,  the  basic  matching  procedure  combines  differences  in  all 
feature  values  which  are  weighted  to  account  for  the  difference  ranges 
of  values  (small  size  differences  (1000  pixels)  are  not  as  important 
for  large  regions  but  small  changes  in  orientation  (0.5  radian)  are 
very  significant) .  Additionally,  differences  in  the  number  of 
relations  in  the  model  and  the  number  of  the  same  relations  between 
the  corresponding  image  segments  are  used.  All  of  these  components 
are  given  a  strength  (high,  medium,  low)  to  control  their  impact  on 
the  final  match  result  (i.e.  features  known  to  be  marginally  useful 
are  given  a  low  strength,  and  those  considered  very  important  are 
given  a  high  strenth) .  In  the  earlier  implementation  the  absolute 
value  of  the  rating  ranged  from  0  upward  to  10000  or  more.  These 
values  are  converted  to  the  range  (0.0,  1.0]  by  using  the  reciprocal 
of  the  value  plus  1. 

The  stochastic  labeling  operation  uses  the  matching  procedure  for 
two  distinct  purposes.  First,  it  is  necessary  to  determine  the 
initial  likelihood  of  a  particular  assignment  (i.e.  a  rating  without 
consideration  of  neighbors,  or  in  other  words  use  only  the  unary 
relations).  Second,  the  compatibility  of  particular  assignments  for 
two  objects  must  be  computed  (i.e.  the  rating  using  the  relations 
between  two  objects) . 


The  computation  of  initial  likelihoods  cannot  use  relations 
between  objects  since  their  use  depends  on  assignments  of  model 
elements  to  image  segments.  Therefore,  the  initial  match  is  limited 
to  feature  values  (color,  size,  shape,  etc.).  A  model  element  is 
compared  with  all  the  image  segments  using  our  basic  matching 
procedure.  The  best  matching  segments  are  kept  for  further  analysis, 
currently  up  to  thirty  are  used  or  up  to  the  point  where  the  worst  is 
1/10  of  the  best  (whichever  given  the  smaller  set).  Then  the  match 
results  are  scaled  so  that  they  sum  to  1,  to  be  treated  as 
probabilities  in  the  stochastic  labeling  procedure. 

Clearly,  if  feature  values,  alone,  are  sufficient  to  locate 
correct  matches,  then  the  process  could  stop  at  this  point.  But,  even 
though  features  are  sufficient  for  some  well  defined  objects,  they  do 
not  locate  most  correspondences  by  themselves. 

After  an  initial  application  of  the  relaxation  procedure  several 
assignments  may  be  made.  At  this  point,  the  computation  of  initial 
likelihoods  for  an  object  can,  and  does,  use  the  relations  with 
assigned  elements.  This  means  that  initial  likelihoods  are  always 
computed  with  all  available  information,  initially  only  feature  values 
then  an  increasing  number  of  relations.  Therefore  successive  steps 
which  are  using  more  information  can  more  reliably  match  the  less  well 
defined  objects. 

The  compatibility  measure  computes  the  effect  of  making  an 
assignment  for  one  element  on  the  assignment  for  another  element.  The 
interaction  between  objects  and  their  assignments  is  through  the 
relations  between  them,  so  that  the  compatibility  measure  is  based  on 
these  relations.  Here  again,  the  same  basic  matching  procedure  is 
used  to  compute  how  well  the  relations  in  the  model  match  the 
relations  in  the  image.  The  procedure  has  been  modified  so  that  only 
the  relations  between  the  two  model  elements  and  between  the  potential 
corresponding  image  elements  are  considered. 


In  some  implementations  of  a  stochastic  labeling  procedure,  all 
the  possible  compatibility  measures  are  computed  once  at  the 
beginning,  but  because  of  the  total  number  of  possibilities  which 
would  be  required,  this  can  not  be  done.  Since  only  a  small  fraction 
of  the  total  number  are  ever  required,  they  are  computed  as  needed. 
Clearly,  in  one  experiment,  the  same  value  will  be  computed  several 
times,  but  the  cost  is  small. 

The  compatibility  measure  is  computed  using  the  most  likely 
assignments  for  the  second  object.  The  individual  matches  are 
weighted  by  the  likelihood  of  a  particular  assignment  so  that  more 
likely  assignments  contribute  more  to  the  result.  The  number  of 
likely  assignments  to  be  considered  is  determined  by  an  input 
parameter.  The  greater  the  number  of  assignments,  the  greater  time 
the  procedure  will  take.  Experiments  have  indicated  that  using  more 
than  one  or  two  assignments  does  not  contribute  much  to  the  final 
results  (see  the  final  section) .  If  an  object  has  a  firm  assignment 
(i.e.  some  segment  has  been  selected  as  corresponding  to  a  model 
object)  then  this  assignment  is  included  in  the  compatibility 
computation  in  addition  to  the  number  of  assignments  specified  by  the 
parameter  described  above. 

Relaxation  Algorithm 

This  section  describes  the  basic  stochastic  labeling  and  the 
optimization  techniques  algorithms.  The  second  part  of  the  section 
presents  the  variations  of  the  basic  procedure  which  were  required  for 
this  symbolic  matching  problem. 

General  Description 

Relaxation  labeling  attempts  to  efficiently  solve  a  very  general 
problem  in  Pattern  Recognition  and  Artificial  Intelligence:  given  a 
set  of  units  U  and  a  set  of  names  N,  assign  names  to  units  given 
actual  measured  features  and  a  world  model .  There  are  two  broad 
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classes  of  Relaxation  techniques.  The  first  one,  called  discrete 
Relaxation,  handles  the  case  where,  for  every  unit,  we  can  know  if  a 
name  is  possible  or  impossible.  Discrete  Relaxation  is  then  an 
efficient  way  to  solve  the  search  problem  of  finding  a  labeling  of  the 
units.  Continuous  Relaxation  handles  the  cases  where  we  know  more 
than  just  whether  names  are  possible  or  impossible,  namely  a  measure 
of  their  likelihood. 


In  the  first  case,  the  world  model  consists  of  a  binary  relation 
RC(UxN) x (UxN)  that  determines  whether  assigning  name  n^  to  unit  u^  is 
compatible  with  assigning  name  n2  to  unit  U2 •  In  the  second  case, 
this  compatibility  is  given  by  a  positive  number  c (u^ ,n^ ,u2 ,n2 ) ,  which 
is  small  if  the  compatibility  is  weak  and  large  if  it  is  strong. 
Function  c  is  defined  in  general  only  over  a  subset  S  of  (UxN)x(UxN). 
To  every  unit  u.  and  name  n,  we  can  thus  associate  the  set  V. (k)  of 
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related  units  u^  such  that  there  exists  a  name  n^  for  which 
(u±  ,1^  ,Uj  ,n£)  is  in  S. 


In  this  paper,  we  are  solely  concerned  with  continuous 

Relaxation,  where  for  every  unit  u. ,  there  is  a  corresponding 

t1 

probability  vector  p^  «  [p^  (1 ) .  ,p^  (L^  )  ]  where  p^(k)  (l_<k£L^)  measures 

the  probability  that  unit  u^  has  the  name  n^ .  The  set  of  all  vectors 
is  called  a  stochastic  labeling  of  the  set  of  units  U.  As  proposed 
in  [7,8]  we  can  also  define  for  each  unit  u^  a  compatibility  or 
prediction  vector  =  [q^  (1 ) . ,  .q^  (L^  )  1  that  tells  us  what  p^  should 
be,  given  the  probability  of  assignments  p^  at  neighboring  units  u^ 
and  the  world  model  embedded  in  function  c.  For  simplicity  we  will 
rewrite  c  (u^,nk,Uj  ,n^)  as  c(i,k,j,Jt)  in  what  follows. 

As  described  in  [7,8]  we  can  take: 


qA  (*> 


CKfk) 


L. 

EQim 

1=1 


(1 ) 


where 
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>iOO  =  IE  c(i,k,  j,A)Pj  (i.) 

Uj  inVi(£)1  in  wj 


and  W.  is  a  subset  of  the  possible  names  for  unit  u..  In  the 
D  3 

experiments  reported  in  [7,8]  we  took  Wj  =  [l,...,Ljb  but  because  of 

the  large  number  of  possible  names  in  this  task  and  for  the  sake  of 

efficiency  we  took 


Wj  =  (set  of  n  most  likely  names} 

usually  with  n  *  1.  That  is,  for  every  neighbor  of  a  unit  we 
considered  only  the  contribution  of  the  most  likely  name  in  Eq.  (2). 

In  [7,8]  we  proposed  to  use  local  measures  of  consistency  and 
ambiguity  of  the  form 


«l  |jj  +  (1-“)*^  (?) 


where  11*112  *s  fche  usual  Euclidean  norm,  the  quadratic  entropy 
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and  a  a  weighting  factor  adjusting  the  relative  importance  of 
consistency  versus  ambiguity.  It  was  found  in  [9]  that  an  even  better 
measure  is  given  by  the  inner  product 


The  global  measure  is  then  an  average  over  the  set  of  units  of  the 
local  measures.  We  can  thus  define 
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all  units 
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as  a  global  criterion  over  the  set  of  units  that  measures  the 
consistency  and  ambiguity  of  the  labeling. 
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The  labeling  problem  is  now  equivalent  to  the  following: 

given  an  initial  labeling  p|°^  ,  find  the  local  maximum  of 
criterion  C  closest  to  the  original  labeling  subject  to  the 
constraint  that  vectors  p^  are  probability  vectors. 


This  is  typically  a  constrained  optimization  problem  which  as  shown  in 
[7,8]  can  be  efficiently  solved  by  using  steepest  descent  techniques. 
In  particular,  we  can  attach  to  every  unit  u.  a  local  gradient  vector 
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and  define  an  iteration  scheme  as: 
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Sfn+1)  =  p.(n)+  d  P(g.(n)  }  n>  0 
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where  Pn  is  a  positive  number  and  P  a  linear  projection  operator.  For 
a  description  of  P  see  [7], 

It  can  be  easily  shown  that: 
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The  first  term  q^(k)  in  Eq.  (6)  corresponds  to  the  simple  maximization 
of  the  product  p^»q^  in  global  criterion  C,  whereas  the  second 
term  corresponds  to  the  coupling  between  units  through  the 
compatibility  relations.  The  algorithm  defined  by  Eq.  (5)  will  allow 
us  to  evolve  from  the  initial  stochastic  labeling  toward  a  less 
ambiguous  and  more  consistent  labeling. 

Least  Commitment  Versus  Speed ,  Multiple  Matches 

The  task  of  matching  a  model  with  a  symbolic  representation  of  an 
image  obtained  by  a  automatic  segmentation  procedure  presents  the 
important  characteristic  that  matches  may  not  be  unique.  For  example, 
if  a  highway  in  the  image  has  been  separated  by  the  segmentation 
procedure  into  several  disconnected  pieces,  then  each  piece  is  a 
potential  correct  match  for  the  node  "highway"  in  the  model. 

One  possible  way  of  handling  this  problem  would  be  to  relax  the 
constraint  that  vectors  are  probability  vectors  and  interpret  them 
as  confidence  vectors  for  which  every  component  may  vary  between  0  and 
1.  This  would  have  the  benefit  of  delaying  any  firm  commitment  as 
much  as  possible,  but  would  also  have  a  general  tendency  to  slow  the 
convergence  of  the  iterative  scheme.  This  is  why  we  looked  for  an 
alternative  between  an  increase  in  speed  (and  therefore  of  the 
probability  of  making  errors)  and  the  principle  of  least  commitment. 

We  therefore  introduced  the  notion  of  a  Macro-iteration  composed 
of  several  of  the  iterations  defined  by  Eq.  (7)  (Micro-iterations! 
after  which  decisions  are  made  to  assign  names  to  units  based  upon  the 
comparison  of  the  components  of  the  vectors  p^  to  a  threshold  (usually 
80%).  If  one  component  is  larger  than  the  threshold  then  it  is  set 
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equal  to  1  and  unit  u^  is  considered  to  be  assigned  the  corresponding 

name  n,  . 
k 

The  process  is  then  reinitialized  by  computing  new  initial 
probabilities  p^^  .  For  units  u.  which  have  been  assigned  names  these 
are  not  actually  probability  vectors  any  more.  The  sum  of  their 
components  adds  up  to  k^+1  where  k^  is  the  number  of  assigned  names. 
Since  relations  are  used  to  compute  initial  probabilities  only  when 
units*  are  assigned  (otherwise  only  features  are  used) ,  this  also  has 
the  advantage  of  improving  the  original  estimates. 

In  the  optimization  problem  (<?)  the  constraints  are  now: 

L. 
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Pi(4)  =  1  for  all  names  £  which  have 

already  been  assigned  to  unit  ui 

That  is,  we  do  not  allow  the  confidence  value  for  already  assigned 
names  to  be  changed  from  1.  f 


Comparison  With  Other  Approaches  And  Experimental  Results 


We  compared  our  approach  with  results  from  two  other  techniques: 
the  relaxation  procedure  originally  introduced  by  Rosenfeld,  Hummel 
and  Zucker  [10]  (algorithm  RHZ)  and  the  sequential  matching  procedure 
developed  by  Nevatia  and  Price  [6]  (algorithm  NP) .  The  nonlinear 
iterative  updating  formula  of  algorithm  RHZ  can  be  expressed  as 


(n+l  pf’oOQ^k) 
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where  the  Qi's  are  computed  in  the  same  way  as  in  Eq.  (2).  On  this 
particular  proble,  algorithm  RHZ  has  a  tendency  to  converge  toward 
ambiguous  solutions.  The  reason  is,  of  course  that  this  algorithm 
does  not  take  into  account  completely  the  coupling  between  units  as 
reflected  by  the  c ( i , k, j ,£) 1 s  and  fails  to  account  for  the  notion  of 
ambiguity.  As  shown  in  Eqs.  (6)  and  (8)  this  is  not  the  case  with  the 
algorithm  described  in  this  paper. 

The  basic  comparison  procedure  used  here  is  the  same  as  used  in 
algorithm  NP,  so  that  many  of  the  results  are  identical.  The  major 
problems  with  the  sequential  method  are  the  lack  of  a  consistent 
method  for  the  handling  of  multiple  assignments  and  errors  which  are 
caused  by  the  order  in  which  the  objects  are  selected  for  matching. 
For  example,  if  a  pair  of  identical  objects  are  near  each  other,  the 
sequential  procedure  can  easily  find  the  wrong  assignment  for  the 
first  object,  thus  forcing  the  second  to  also  be  incorrect,  due  to  the 
limit  of  assigning  one  name  to  only  one  unit  (units  may  be  assigned 
several  names).  Since  the  relaxation  procedure  is  considering  pairs 
of  objects,  the  correct  assignments  advance  to  the  top,  if  there  are 
relations  connecting  the  two  objects. 

Results 


We  have  applied  this  procedure  to  several  different  scenes  or 
portions  of  scenes  and  will  present  results  from  two  of  these  scenes. 
The  first  is  a  high  altitude  aerial  image  with  a  few  major  regions  (a 
city  and  rural  areas)  and  a  number  of  linear  features  (a  major 
highway,  river  channels,  and  roads  along  the  channels),  see  Fig.  1. 
The  second  is  a  closer  view  (of  a  different  area)  ,  with  2  roads 
included  in  the  model.  The  roads  are  described  in  sections  because  of 
the  tendency  of  the  linear  feature  extraction  procedure  to  break  long 
linear  features  into  shorter  segments,  see  Fig.  2. 
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To  illustrate  how  the  relaxation  procedure  operates.  Fig.  3  is  a 
graph  of  the  probabilities  of  various  assignments  for  one  object 
through  a  series  of  macro  iterations. 

Figs,  4  and  5  show  another  aspect  of  the  operation  of  the 
relaxation  process.  The  15  most  likely  assignments  for  one  model 
object  are  shown  with  the  most  likely  one  labeled  1  and  least,  likely 
labeled  15.  The  first  picture  is  the  initial  likelihoods  based 
primarily  on  feature  values,  the  succeeding  pictures  are  the  likely 
assignments  after  2  and  5  (micro)  iterations.  Figs,  f  and  7  show  the 
final  results,  with  the  objects  outlined  and  labeled.  Many  of  the 
labels  overlap  since  adjacent  linear  features  are  being  presented. 
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is  labeled  " SOUTH-HIGHWAY "  in  the  final  results  of  Fig.  6.  The 
tick  marks  along  the  horizontal  axis  indicate  each  Macro  iterat 
The  three  segments  which  rise  to  1.0  are  all  correct  assignment 
Note  also  that  the  initial  best  assignments  is  not  considered 
after  the  firm  assignment  made  on  the  second  macro  iteration. 


Sequence  of  15  most  likely  labels  for 
a  river  area.  a)  Ordered  by  initial 
probabilities.  The  correct  assignment 
is  also  the  >most  likely,  b)  Likely 
assignments  after  2  iterations. c)  Likely 
assignments  after  4  iterations  (now  only 
10  with  nonzero  likelihood) . 
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(c) 

Fig.  5  Sequence  of  15  most  likely  labels  for  a 
road  segment.  a)  Ordered  by  initial  pro¬ 
babilities  based  primarily  on  feature 
values.  The  correct  label  is  the  sixth  most 
likely  label.  b)  Likely  assignments  after 
3  iterations,  c)  Likely  assignments  after  6 
iterations . 
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2.2  Decomposition  and  Decentral i zation  Techniques  in  Relaxation 
Label ing 

O.D.  Faugeras 


Introduction 


One  of  the  most  important  tasks  in  Image  Analysis  is  to  assign 
names  to  objects  present  in  the  scene.  This  semantic  analysis  usually 
happens  after  some  symbolic  analysis  like  segmentation  and  shape 
analysis  has  been  performed  but  it.  does  not  have  to.  Low  level  vision 
can  also  be  considered  as  semantic  analysis  when  for  example  we 
attempt  to  decide  whether  an  edge  is  present  or  not  at  some  point,  in 
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an  image  (names  are  here  edge/noedge)  or  when  we  try  to  segment  a 
picture  into  objects  and  background  (names  are  here 
objects/background) .  Usually  the  task  of  assigning  names  to  objects 
is  very  difficult  on  the  basis  of  only  features  measured  for  that 
object.  This  is  typically  a  Pattern  Recognition  problem  and  errors 
are  likely  to  be  made  for  reasons  such  as  noisy  data  and  inadequate 
feature  extraction.  One  of  the  key  ideas  in  getting  around  that 
problem  is  to  delay  any  firm  commitment  until  all  the  available 
contextual  information  has  been  used. 

Waltz  [1]  was  one  of  the  first  to  use  the  notion  of  constraints 
induced  by  a  world  model  to  solve  the  problem  of  analyzing  a  scene 
made  up  of  complex  polyhedra.  His  ideas  were  generalized  by  Rosenfeld 
et  al.  (2]  who  proposed  a  completely  parallel  version  of  Waltz 
labeling  algorithm.  Davis  and  Rosenfeld  presented  some  appl ications 
of  the  so-called  "discrete"  relaxation  techniques  to  template  matching 
[3],  shape  matching  [4]  and  waveform  parsing  T5],  Haralick  r*l 
introduced  a  general  framework,  the  theory  of  arrangements,  in  which 
the  solution  to  some  scene  analysis  tasks  can  be  found  by  constructing 
a  homomorphism  from  one  arrangement  to  another.  To  help  determine  all 
these  homomorphisms,  he  generalizes  to  N  dimensions  the  Waltz 
filtering,  the  discrete  scene  labeling  relaxation  and  the  network 
consistency  relation  r 7 ] . 

Barrow  and  Tenenbaum  f8]  and  Rosenfeld,  Hummel  and  Zucker  r?l 
introduced  the  idea  of  probabilistic  or  stochastic  labeling  where  we 
not  only  know  possible  and  impossible  classes  for  every  object  but 
also  a  measure  of  their  likelihood.  Similarly  the  world  model  is  now 
described  continuously  in  terms  of  compatibility  coefficients  or 
conditional  probabilities  between  labels.  The  nonlinear  algorithm 
proposed  in  (2]  has  been  used  in  different  applications,  edge  and  line 
enhancement  [9,10,11],  clustering  Til, 12],  stereo  vision  and  movement 
detection  [13],  Similar  ideas  have  been  used  by  Marr  and  Poggio  fl4]. 
Theoretical  analysis  of  the  convergence  and  stability  properties  of 
this  algorithm  has  proven  to  be  difficult  as  shown  by  Zucker 
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*15,16,17].  Faugeras  and  Berthod  [18,19,20]  have  proposed  a  different 
approach  based  upon  explicit  use  of  compatibility  and  ambiguity  to 
define  a  global  criterion  on  the  set  of  objects.  They  showed  that  the 
criterion  could  be  minimized  using  projected  gradient  techniques.  A 
similar  idea  has  been  proposed  by  Ullman  r 2 1 ] .  We  will  not  proceed  to 
describe  further  the  problem  of  stochastic  labeling  viewed  as  an 
optimization  task . 

Stochastic  label inq  as  a  problem  i n  Opt imi zat i on 
Definition  of  criterion: 


Formally,  we  are  given  a  set  of  N  objects  a. ,a„ , . . . ,a„  which  fall 

1  Z  N 

into  L  classes  A,,...,AT.  To  every  object  a.  we  assign  a  probability 

i.  Li  1 

noted  p .  ( A,  )  to  belong  to  class  A..  This  is  conveniently  represented 

as  a  vector  p.=  [p. ( A. ) , . . . ,p . ( A,  ) ]  .  Objects  are  also  related  to 
1  1  -L  X 

one-another:  we  will  denote  by  V.(A,)  the  set  of  objects  related  to 

1  K 

object  a.  for  label  \  .  This  means  that  for  every  object  a.  in  V.  (a,  ) 

1  K  J  1  K 

we  are  given  the  conditional  probabilities  p. . (A.  |A.)  that  object  a. 

13  k  y,  1 

belongs  to  class  \  given  that  object  a.  belongs  to  class  A0, 

£=1,...,L.  At  this  point  we  would  like  to  stress  the  fact  that  the 

vectors  p.  are  the  results  of  measurements  taken  on  the  data  to  be 
1 

analyzed  and  that  the  p. .(A,  |A0)'s  constitute  our  a  priori  contextual 

1  ]  ,  K  K  ^ 

information  or  world  model.  Due  to  noisy  data,  vectors  p.  are  not 
always  compatible  with  respect  to  the  contextual  information  and 
ambiguous,  that  is,  in  general  we  do  not  have 


pi<V 


(1  ) 


and  probability  vectors  are  not  unit  vectors  ro, . . .  0 , 1 ,  P. . .  P 1T. 
From  these  two  facts  stems  the  need  for  designing  algorithms  that  use 
contextual  information  and  the  notion  of  ambiguity  to  modify  the 
probability  for  every  object  in  such  a  way  that  consistency  is 
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increased  and  ambiguity  decreased. 


The  approach  that  we  suggested  (18,19,20]  was  based  upon  the 

minimization  of  a  criterion  defined  over  the  set  of  objects.  The 

class  of  criteria  that  we  proposed  was  meant  to  be  a  global  measure  of 

the  amount  of  consistency  and  ambiguity  present  in  the  stochastic 

labeling.  Local  consistency  is  measured  as  the  norm  of  the  difference 

of  the  probability  vector  p.  for  object  a-  and  a  consistency  vector 
->•  rp  1  X 

q-fqiUi) . qiUL>]  whose  coordinates  are  some  normalized  averages 

of  what  we  would  expect  the  probabilities  of  labels  for  object  a^  to 
be  when  taking  into  account  both  the  measurements  at  objects  related 
to  a^  and  the  contextual  information.  More  precisely,  we  have: 


qi(Xk>  “  -E 


E  Qi<M 

4=1  1  x 


m 


and 


W  - 


_  / _ i 

'  (ff5 


FFT 


E 


s*.yy 

i]k  I 


0  <  r  <  °° 


(?) 


The  positive  numbers  s.  ..  are  given  by: 

1  JK 


i-i 

sijk  “  “ijk  £ 


(4) 


The  a. .Js  are  positive  and  weight  the  importance  we  attribute  to  the 

i  J  K 

different  objects  a^  in  V^(Ak)  in  evaluating  the  consistency  of  label 
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A ^  for  object  (index  j)  as  well  as  the  importance  we  attribute  to 
that  particular  label  in  ascertaining  the  total  consistency  for  object 
ai  (index  k) .  As  was  previously  mentioned,  this  quantity  is  measured 
by  the  vector  norm 

L  i. 

llPi-Ma"  (£  lplUk>-^i(xk>l“)“  0  <  a  <  - 


Global  consistency  is  then  measured  by  an  average  taken  over  the  set 
of  objects  of  the  local  consistency  measures,  for  example  the 
arithmetic  mean: 


C 


N 


k  S  IlMil'a 

x=l 


(«) 


In  the  same  spirit,  local  ambiguity  can  be  measured  riR, 19,201  by  the 
entropy  of  the  probability  vector  p^  as  defined  by 

L 

Hi  ”  2  pi(xk)d-pi(xkn  =  i-  tiPilli  e) 

Global  ambiguity  is  then  measured  by  the  average  of  the  local 
ambiguity  measures: 

N 

A  =  1-jj  II  P±  II  2  (*> 

The  total  criterion  J  is  then  defined  as 


J  =  cjC+  (l-w)A 


(9 ) 


where  coefficient  w  which  varies  between  0  and  1  weights  the  relative 
importance  we  attribute  to  consistency  versus  ambiguity. 


Computational  aspects 

Let  us  denote  by  v  the  vector  of  RM=RLx . . . xRL  (M=NL)  equal  to 

N  times 

(Pl»i?2  r  •  •  •  »Pjg)  •  We  can  define  the  stochastic  labeling  problem  as 
follows : 

given  an  initial  stochastic  labeling  v°  find  the 
stochastic  labeling  u  that  corresponds  to  the 
local  minimum  of  criterion  J  (v)  which  is  closest 
to  v°  subject  to  the  constraint  that  if 
u= (p^ , . . . ,pN)  then  p^  is  a  probability  vector  for 
'  i  =  l , . . . ,  N 

The  fact  that  we  are  not  requiring  to  find  the  absolute  minimum  of  J 
corresponds  to  the  desire  to  evolve  toward  a  stochastic  labeling  which 
depends  on  the  initial  measurements.  The  contrary  would  imply  that 
the  final  labeling  would  be  independent  of  the  original  set  of  objects 
and  would  depend  solely  on  the  world  model! 

The  optimization  problem  may  in  practice  be  of  a  very  large 
dimensionality  if  the  number  of  objects  and  possible  labels  is  not 
kept  within  reasonable  bounds.  Unfortunately,  for  many  practical 
problems  this  is  not  possible  and  we  have  to  find  means  of  getting 
around  the  "curse  of  dimensionality".  In  r 2 0 1  we  show  that  a  large 
amount  of  parallelism  can  be  introduced  in  the  minimization  process: 
if  we  attach  to  every  object  a^  a  processor  r^  connected  only  to 
processors  r^  attached  to  objects  a^  related  to  a^ ,  then  the  global 
criterion  J  can  be  minimized  by  having  processors  r^  perform  simple 
operations  mostly  in  parallel  while  a  simple  sequential  communication 
process  allows  them  to  work  toward  the  final  goal  in  a  coordinated 
fashion . 


Another  set  of  related  ideas  can  be  found  in  the  so-called 
Decomposition  and  Decentralization  techniques  which  have  been 
developed  to  solve  similar  problems  in  Economy,  Numerical  Analysis, 
Systems  Theory  and  Optimal  Control  [22,23,24,25].  Decomposition 
techniques  proceed  from  an  algorithmic  standpoint:  we  are  confronted 
with  a  problem  of  large  dimensionality  and  try  to  substitute  for  it  a 
sequence  of  problems  of  smaller  dimensionalities.  Decentralization 
techniques  take  a  different  viewpoint:  we  are  confronted  with  a  global 
problem  and  have  at  our  disposal  P  decision  centers.  The  question 
that  we  ask  ourselves  is  whether  it  is  possible  to  solve  the  global 
problem  while  letting  the  decision  centers  solve  only  local  problems. 

Our  purpose  here  is  to  show  that  the  structure  of  our  criterion 
allows  us  to  develop  both  types  of  techniques.  In  order  to  do  that  we 
need  to  develop  a  few  notations. 

Grouping  objects  into  sets 

Let  k=l,...,P  be  nonoverlapping  sets  of  objects,  that  is 

nk°  =  0  if  k  /  l  (10) 

For  ever  object  a^ ,  i=l,...,N  we  will  denote  by  \L  the  set  of  all 

objects  a.  related  to  it  for  some  labels  : 

J  K 


v. 


L 

U 

k=l 


w 


(11) 


Two  sets  and  can  "interact"  if  for  some  object  a^  in  ,  some  of 
its  neighbors  are  in  that  is  if  VDfi^/o.  We  will  then  define^  to 
be 
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p 

=  u 
4*1 


1?  k 


all  objects  in  fto  which  are  neighbors)  Uy] 
of  some  objects  m  ft,  I  k 


and  N^  to  be  the  number  of  objects  in  ft^.  Notice  that  we  may  now  have 
for  some  pairs  (k,  i)  with  kfl.  Let  us  cal  ] 


min  (k,  4) max  (k,  l)  r  *'min  (kf  i)max  (k,  Jl) 

the  number  of  objects  in  it.  Given  a  subset  of  objects  we  will 

denote  by  vL  the  restriction  of  the  vector  v=  (p-^ , . . .  ,pN)  to  objects 

in  k 


the  previous  intersection  and  N 


Let  us  now  take  a  closer  look  at  the  structure  of  criterion  J. 
Clearly  it  can  be  rewritten  as  a  sum: 


N 


J(v)  =  £  JJv) 
k=l  K 


(13) 


where,  according  to  equations  (6),  (8)  and  (9) 


Jk(v)  =  ^  <w||pk-qklla+ (i-m)  (i -IIp^II^  > 


(14) 


From  the  definition  of  vector  qk ,  it  is  seen  that  criterion  Jk  (v) 
depends  solely  upon  vectors  Pj  for  objects  aj  in  vk  or  to  put  it 

the  restriction  of  v  to 


is 


another  way  that  J  (v)=J  (v.)  where  v 
—  (0 )  (u )  k  k  K 

objects  in  (ft^  ={a^}.  This  decomposition  of  criterion  J  as  a  sum 

of  N  criteria  corresponds  to  the  partition  of  the  set  of  objects  into 

the  subsets  ft^?  ,  k=l,...,N.  More  generally,  to  every  partition  ft  , 

A  K 

k=l,...,P  of  the  set  of  objects  we  can  associate  a  decomposition  of 
criterion  J  as  a  sum  of  P  criteria: 
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where 


J  (v)  =  Ec,*) 
k=l  K 


Ck(v) 


E 

i  in  ft, 
k 


J.(V)  S  Ck(vk) 


(16) 


and 


Decomposition  techniques 

Our  problem  is  to  minimize  a  criterion  J  (v)  given  by: 

P  P 

J(v)  =  23cv(v)  E  23c,(vw) 
k=l  K  k=lK 


(17) 


(IB) 


where  v=(j$  )  is  subject  to  the  constraint  that 

1  N 

probability  vector  for  i=l,...,N  and 


Ck(v) 


(1?) 


The  set  of  constraints  can  be  expressed  as 
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■i 


TH 


(20) 


P 

k  =  n 

k=l 

Kk={  v|  v=  (p1# . . .  ,pN)  such  that  for  all  objects  a^  in  Q^p^  is  a 
probability  vector). 

The  decomposition  of  K  according  to  (20)  corresponds  to  a  splitting  of 
constraints:  the  idea  is  to  minimize  a  criterion  with  a  lesser  number 

of  constraints. 


Sequential  algorithm 


Let  Tl' *  *  *  3  sec?uence  ^  of  strictly  positive  numbers  we 

define  a  family  of  vectors  un+P  i=l,. 

any  u ^  in  K  and  assuming  that  we  have 
>  n+k 

u  F  is  the  element  of  K k  that  minimizes 


.  .  , P,  n=0 ,  1 ,  .  .  . 
computed  u  . . 


starting  from 

^n-4-k-  1 

then 


-+n+__ 

,u  , 


Ck<V>+2T 


n+ 


k-1 


v-u 


(21 ) 


We  can  associate  with  the  sequence  of  vectors  tP+p  the  sequence  u>N+p 
given  by 


-N+| 


U) 


■’a 


N  ,k 

1  V  -*n+P 
=  —  2_J  T  U 

0N  n=0  n 


(22) 


where  n  is  the  sum  of  the  numbers  t  ,  n=0,...,N. 

N  n  k 

It  can  be  shown  [25]  that  the  limit  when  N  goes  to  infinity  of  u  pis 
a  local  minimum  of  criterion  J  for  k=l,...,P.  The  algorithm  then  goes 
as  follows: 


Algorithm  DS : 

step  1.  (Initialization)  start  with  u®  in  K,  set  n  equal  to  0 
step  2.  (Iteration  initialization)  set  k  equal  to  1 
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closest  to 


step  3. 


un  P 


problem)  find 
of  criterion 


the 


minimum 


in  Kk 


Ck(v)+TT 


II  v-u 


.  n+ 


k-1 


n 


(23) 


step  4. 

(Test 

equal 

for  end  of  loop) 

P  go  to  step  3. 

set  k  equal  to 

k+1  . 

If  k 

less  or 

step  5. 

(Test 

equal 

for  convergence) 

to  n+1.  and  go  to 

if  converged 

step  2. 

stop 

else 

set  n 

If  we  stopped  for  n=N  the  result  is 


(?a ) 


for  k=l,...,P.  In  practice  of  course,  since  we  stop  after  a  finite 

^ 

number  of  iterations,  vectors  J  p  are  likely  to  be  different  for 
k=l,...,P.  They  should  be  nonetheless  close  to  one  another  and  we  can 
either  choose  one  of  them  as  the  final  result  or  their  average 


r?F) 


Parallel  a  Igor i thm 


a 


0 


Again  we  start  with  a  vector  u^ 
n  i-  k 

..,u  ,  we  define  u  P+1  for  k  =  l  , 


in  K  and  supposing  that  we 
..,P  as  a  local  minimum  of 


know 
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(26) 


CR(v)  =  2T-||^n|P 
n 


for  v  in  Kk*un+^  is  then  defined  as 


un+1  -  j 

k=l 


(27) 


Just  as  for  the  sequential  algorithm,  it  is  possible  to  show  [25] that 
the  sequence  of  vectors  u>  P+I  defined  as 


k  N  ^  k 

,N+P+I  _  1  V  -n+P+T 

J  a.  Tnu 

N  n=0 


(28) 


has  a  limit  when  N  goes  to  infinity  and  that  this  limit  is  a  local 
minimum  of  criterion  J  for  k  =  1,...,P.  The  algorithm  then  goes  as 
follows : 


Algorithm  DP: 

step  1.  (Initialization)  start  with  u° 
step  2.  (Solve  problems  in  parallel) 
closest  to  un  of  criterion 


in  K  and  set  n  equal  to 

.  .  -*tih — — — 

find  the  minimum  u  P+1 


0 

in 


Ck(v)+2T 


i '+  ->ni 

v-u 


n 


(29) 


for  k  =  1, . . . ,P 


Step  3.  (Test  for  convergence)  if  converged  stop  else  set  un+^ 

1  ^  ->-n+P+T 

equal  to  p  Z  u  ,  n  equal  to  n+1  and  go  to  step  2, 

F  k=l 
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If  we  stopped  at  N  the  result  is 


CO 


-*n+ 

u 


k 

P+1 


n 


(30) 


for  k  =  1,...,P.  Just  as  in  the  case  of  the  sequential  algorithm  since 

•  •  — 
we  stop  after  a  finite  number  of  iterations,  vectors  to  P+1  are  likely 

to  be  different  for  k  =  1,...,P.  The  same  remedy  can  be  used. 

Simple  case  analysis 


Let  us  consider  a  set  { a^ ,a2 ,a^ , a4 )  of  four  objects  related  as 

in  figure  1,  that  is  to  say 


aQ  =  a^  and 

VL  =  i  =  1,...,4  where  '  4 

a5  =  ax 

Let  us  now  choose  for  example 


«i  =  { a i >  i  =  1, . .  .  ,4 

We  thus  have  TL  =  {a^}  U  for  i  =  1,...,4  and 

=  (P3.P2.P4) 

v2  =  (P3.P2.P3) 

V3  "  (P2.P3.P4) 
v4  =  (Pi.P3.P4) 

We  also  have: 

J 


(31) 


V?)  =  CjfJj)  =  fllJj-Sjl2*  i^d-llpjl2) 

where  according  to  figure  1  and  equations  (2-4),  vector  <5^  depends 

solely  upon  vectors  py  and  p . .  Starting  with  initial  conditions 

-+0  .  -*0  ->0  0  -+0  %  • 

u=(P^#p2»P3»  P4) t  in  the  case  of  the  sequential  algorithm  we  minimize 
first 

cl(J)+jl_(||fl-fJ||2+||g2-p“||2+||p3-p°|l2+l54-5j||2)  ,32) 


since  C^(v)  does  not  depend  upon  p  the  result  will  be 


u°+*  =  <2*  &  -°  -i 

(Pl r P2  > P3  »  P4 ) 


(3?) 


Similarly  we  will  have: 


0+ 


2  2  2  1 


UT-r  *r-  t  t  T 

-*■  4  ,-*A  -+4  -+4  ->4  . 

u  —  (P^  »?2  »P3 'P4  ^ 


(34) 


+0+J  +0+* 

In  the  case  of  the  parallel  algorithm  we  first  compute  u  =u  but 

we  then  have 


0+2  222 
-*-0+B'  -*.5  -+5  -*-5  ->° 

u  =  (P^  »P2 'P3 'P4  ^ 

-0+l 

which,  in  general  is  different  from  u 


(35) 


Thus  the  two  algor ithms  are  not  trivially  equivalent  in  the  sense 
that  at  every  iteration  they  compute  essentially  the  same  thing  in 
different  ways.  They  are  nonetheless  equivalent  in  a  deeper  sense 
that  in  the  limit  they  converge  toward  the  same  answer.  In  practice 
it  would  be  interesting  to  know  whether  or  not  one  algorithm  converges 
faster  than  the  other. 


We  would  like  to  point  out  that  this  simple  example  is  only  meant 
to  help  the  reader  get  a  more  concrete  feeling  for  the  concepts 
developed  in  the  previous  sections.  Because  of  the  small  number  of 


37 


objects  in  the  initial  problem  we  replaced  a  problem  of  dimension  4l 
with  4  problems  of  dimension  3L.  Of  course  the  gain  in  dimensionality 
can  be  a  lot  more  for  larger  size  problems  as  encountered  in  practice. 
The  question  of  how  to  optimally  choose  numbers  xn 1 s  to  speed  up 
convergence  is  largely  unsolved  although  some  preliminary  results  have 
been  obtained  [26]. 

Price  decentralization  techniques 


General  Algorithm: 

The  idea  here  is  to  let  the  decision  center  solve  local  problems 
and  coordinate  them  in  such  a  way  that  the  overall  result  is  the 
minimization  of  the  global  criterion  J.  When  minimizing  C  (v.  )  over 
ft^  the  k  decision  center  is  going  to  compute  vectors  for  objects 

a  •  in  regions  ft.  adjacent  to  region  ft,  and  there  is  no  reason  for 
having  Pfcj=Pjij'  P^j  being  the  result  of  the  minimization  of  C^fv^)  by 
the  decision  center  over  .  Thus  the  need  for  coordination  which 
is  achieved  by  attributing  a  price  to  the  discrepancy  between  p,  .  and 

p£j  * 

Formally,  we  need  to  consider  intersections  between  sets 

ft^  and  ft^  (k<£)  and  define 


H' 


P-1 

n 

k»l 


P 

n 


(RL)Nkil 


JL=k+l 


(3M 


~  L  k 

For  each  set  ft^  we  then  define  a  linear  application  Bk  of  (R  )  into 

4*  as : 
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sk'vk; 


where,  as  usual,  v. 


I. 


S  i 

The  global  problem  J 


|* 

'Lkj 

j  >  k 

A  I* 

|2-ik 

i  <  k 

(37) 

0 

otherwise 

is  the  restriction  of  v^  to  objects  in 

min  J (v) 

(3R  ) 

is  now  equivalent  to  the  problem: 


mm 


P 

E 


V -j  f  •  •  •  *  V 


*  k=l 


ck(V 


(  subject  to  XlBk(vk)  =  <5 


k=l 


P9) 


which  in  turn  is  equivalent  to  finding  a  saddle  point  of  the 
functional  (Lagrangian) : 


L<v, 


» v_ , r ) 


(40) 
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where  the  price  vector  (Lagrange  multiplier)  r  is  a  vector  of  y  and 
denotes  the  inner  product. 

The  algorithm  then  goes  as  follows: 

Algorithm  PD: 

step  1.  (Initialization)  set  £  equal  to  1,  choose  a  price  vector 

?(1) 

step  2.  (Local  problems)  solve  the  P  "local"  problems 

Min  Ck(vk)+?U)  *Bk(vk)  (41) 

vk 

the  results  are  k=l,...,P. 

step  3.  (Test  for  convergence)  if  converged  stop  else  set  £ 

equal  to  £+1  and  go  to  step  4. 

-*-(?  +  l) 

step  4.  (Update  the  price)  r  is  obtained  by 

r(l+l)  -  r° l)+PiY,Bk(vll))  (42) 

1 

where  is  an  acceleration  parameter,  go  to  step  2. 

Just  as  in  the  case  of  the  Decomposition  techniques  we  have  to 
minimize  a  sequence  of  slightly  perturbed  versions  of  criteria  Ck(vk) 
for  k=l,...,P.  Two  important  differences  must  be  pointed  out:  in  the 
case  of  the  price  Decentralization  technique,  the  perturbation  term 
depends  solely  upon  v^,  rather  than  v  for  Decomposition  techniques  and 
is  a  linear  function  of  v^  rather  than  a  quadratic  one. 
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Simple  case  analysis 

Just  as  in  the  case  of  Decomposition  Techniques,  it  may  he 
helpful  at  this  point  to  examine  more  closely  a  concrete  example.  We 
chose  the  same  as  in  a  previous  section,  we  will  need  this  line  to 
distinguish  between  the  components  of  the  different  vectors  vk=v  ^  , 
for  k=l,2,3,4.  We  will  therefore  write:  ^ 


II 

(pll'P12'P14 

-V 

v2  = 

(P21'P22'P23 

-> 

V3 

{p32'P33'P34 

li 

+  > 

(p31,p33'p34 

We  also  have: 


Cl(?)  ,  Cl(V  =  |||p11-qu||2+iiiH).U.||  ?ii||  2,  («, 


where,  according  to  figure  1  and  equation  (24),  vector  qi;Ldepends  only 


on  vectors  p^  and  p  .  We  must  also  consider 


]£.  .  =  fiin  1,3  =  1,2, 3, 4 

rj  J 


1  <  3 


(4/1) 


According  to  a  previous  section,  we  have: 
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{al'a2} 


{  a2 '  a3  ) 


/.  =  {a,, a.)  V*  -  {a,, a-} 

13  14  1  * 

^  ]  ~  {3^  ,a,}  .  ~  .{soia^ } 

^24  1  J  34 


According  to  definition  (37)  we  have 

B^(v^)  =  {p^i • P^2 /P12 'p14 ' P11 ,p14 ' B B ® ® 

B2(v2^  =  ^_P21,-P22'0'0,0'0'P22'P23'P21,P23,0,0) 
b3(v3)  =  (0,0,  “P32  •  "*p34 ' B _p32 ,-p33 ' B  /  0  ,p33  ,p34  ^ 

B4 {^4  >  =  {0,0, 0,0, ”P4i' -p44 '0*0, -p41 ' -p43 ' ~P4  3 ' ~P44  ^ 


4 

So  that,  indeed,  the  condition  B,  (v,  )=0  is  equivalent  to 

k=l  K 


(45) 

(45) 

(47) 


(48) 
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P11  =  P21  =  P31  =  P41 

P12  =  P22  =  p32  =  p42 

P13  =  P23  =  P33  =  P43 

P14  =  P24  =  P34  =  P44 

-*  (  l  ) 

At  iteration  £  of  Algorithm  PD,  let  us  call  r  the  price  vector  and 
v!Z  ,  k=l,2,3,4,  the  results  of  the  (£-l)th  iteration.  The  "iocal" 

problems  are  then: 

!->-  ( 2, ) 

find  the  local  minimum  of  criterion 

closest  to  v(£  ^  subject  to  the  constraint  that  the  p^ ' s 
( i=k-l , k , k+1 )  are  probability  vectors 

Just  as  in  the  case  of  the  Decomposition  techniques,  this  example 
only  intends  to  help  the  reader  get  a  more  concrete  feeling  for  the 
concepts  developed  in  the  previous  section. 

Conclusions 


We  have  presented  in  this  paper  two  methods  of  approaching  the 
problem  of  reducing  the  dimensionality  of  the  Relaxation  labeling 
problem  viewed  as  an  Optimization  task.  The  first  approach  takes  an 
algorithmic  standpoint:  it  replaces  a  large  dimensionality  problem  by 
a  sequence  of  smaller  dimensionality  problems  threaded  by  a  sort  of 
rubber  band  which  forces  each  of  them  to  take  into  account  the  results 
obtained  by  the  others.  The  coordination  problem  is  tackled  directly 
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in  the  second  approach  in  which  a  price  is  dynamically  assigned  to 
eventual  discrepancies  between  solutions  to  the  different  "small" 
problems.  The  problem  of  comparing  the  two  approaches  is  open  but  it 
is  hoped  that  the  rapid  development  of  relaxation  labeling 
applications  involving  large  numbers  of  objects  will  provide  answers 
to  that  question. 
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2.3  Extraction  of  Texture  Primitives 

F.  Vilnrotter*,  R.  Nevatia,  K.  Price 


Introduction 


In  previous  reports,  we  have  described  a  program  used  to  generate 
descriptions  of  natural  textures  11-21.  (These  references  also 
contain  a  discussion  of  other  work  and  a  list  of  related  references!. 
The  edge  and  direction  images  corresponding  to  the  original  texture 
image  serve  as  inputs  to  the  program. 

In  the  first  part  of  the  program  edge  repetition  arrays  are 
produced.  (An  edge  repetition  array  is  in  essence  the  binary  case  of 
a  gray  level  co-occurrence  matrix).  These  arrays  are  calculated  for  6 
directions  (0°,  30",  90^,  120°,  and  lfp1),  for  both  dark  and 
light  intensity  objects,  at  distances  within  a  range  specified  by  the 
user.  A  comprehensive  discussion  of  edge  repetition  arrays  is  given 
in  ( 1  ]  . 

In  the  second  part  of  the  program  the  edge  repetition  arrays  are 
analyzed  to  determine  whether  there  are  predominant  element  sizes  in 


♦Felicia  Vilnrotter  is  supported  by  a  Hughes  Aircraft  Company  Poctorai 
f el lowsh i p . 


48 


any  of  the  6  scan  directions  and  if  so  whether  these  elements  occur  at 
regular  intervals  within  the  image.  The  details  of  this  analysis  are 
presented  in  [2]  and  will  not  be  repeated  here. 

Figure  2  is  the  description  generated  for  the  raffia  (woven  palm! 
texture  sample  given  in  Fig.  1.  The  information  comprising  this 
description  is  1 -d imens iona 1  in  nature.  The  second  dimension  of  a 
primitive  cannot  be  simply  obtained  by  examining  the  like  intensity 
element  sizes  for  other  scan  directions.  For  example,  looking  at 
Fig.  2  one  might  be  tempted  to  say  that  raffia  is  made  up,  in  part,  of 
light  elements  whose  dimensions  are  3  pixels  along  its  horizontal  axis 
and  5  pixels  along  its  vertical  axis.  This  is,  in  fact,  not  the  case. 
Figure  3  is  an  abstract  representation  of  the  basic  taffia  primitive 
cluster.  The  light  element  found  in  the  vertical  scan  direction  is 
the  medium  intensity  block  A  of  Fig.  3,  while  the  light  element  found 
by  the  horizontal  scan  is  the  high  intensity  block  C  in  Fig.  3.  If 
both  dimensions  of  block  A  were  found  then  the  vertical  dimension 
would  be  character i zed  as  belonging  to  a  light  element  while  the 
horizontal  dimension  would  be  described  as  belonging  to  a  dark 
element.  Clearly,  a  two  dimensional  description  of  the  texture  cannot 
be  constructed  by  combining  its  one  dimensional  parts.  However,  the 
computed  description  can  be  used  as  a  starting  point  for  a  two 
dimensional  primitive  search.  Some  other  primitive  extraction 
techniques  are  discussed  in  T  3 ] . 

Primitive  Extraction  Process 

Motivation  -  There  are  three  strong  indications  of  element  size 
provided  in  the  texture  description  given  in  Fig.  2.  Each  of  these 
predominant  element  sizes  is  the  location  of  a  strong  peak  in  an 
element  size  edge  repetition  array  for  the  raffia  image  in  Fig.  1. 
Knowing  the  exact  locations,  within  the  original  texture  image,  where 
the  edge  matches  contributing  to  these  strong  peaks  occur  would  be 
useful.  It  would  then  be  possible  to  isolate  the  uniform  intensity 
regions  or  textural  elements  being  measured.  Analysis  of  the  set  of 
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Figure  la.  Raffia  Subwindow 


NUMBERS  APPEARING  IN  PARENTHESES  ARE  SCALE  DEPENDENT 
FILENAME  =  RAFFIA -NR1D 
DARK  OBJECT  DESCRIPTIONS 
HORIZONTAL  SCAN  DIRECTION 

NO  EVIDENCE  OF  PERIODICITY  OR  PREDOMINANT  ELEMENT  SIZE 
3D  DEGREE  SCAN  DIRECTION 

NO  EVIDENCE  OF  PERIODICITY 

UEAK  EVIDENCE  OF  PREDOMINANT  ELEMENT  SIZE  (25. MQ) 

VERTICAL  SCAN  DIRECTION 

THERE  IS  STRONG  EVIDENCE  OF  PERIODICITY  (ELEMENT  SPACING  fi-OO) 

THERE  IS  STRONG  EVIDENCE  OF  PREDOMINANT  ELEMENT  SIZE  (3. DO) 

WITH  MODERATE  SUPPORT  FOR  ELEMENT  SPACING  (6- DU) 

RATIO  OF  SIZE  TO  PERIOD  IS  -3fl 

LIGHT  OBJECT  DESCRIPTIONS 
HORIZONTAL  SCAN  DIRECTION 

NO  EVIDENCE  OF  PERIODICITY 

STRONG  EVIDENCE  OF  PREDOMINANT  ELEMENT  SIZE  (3-00) 

3D  DEGREE  SCAN  DIRECTION 

NO  EVIDENCE  OF  PERIODICITY  OR  PREDOMINANT  ELEMENT  SIZE 
VERTICAL  SCAN  DIRECTION 

THERE  IS  STRONG  EVIDENCE  OF  PERIODICITY  (ELEMENT  SPACING  fi-OO) 

THERE  IS  STRONG  EVIDENCE  OF  PREDOMINANT  ELEMENT  SIZE  (5. DO) 

WITH  MODERATE  SUPPORT  FOR  ELEMENT  SPACING  (fi-DD) 

RATIO  OF  SIZE  TO  PERIOD  IS  -L3 

bO-,120  AND  150  DEGREE  SCAN  DIRECTIONS  FOR  BOTH  LIGHT  AND  DARK  OBJECTS 
NO  EVIDENCE  OF  PERIODICITY  OR  PREDOMINANT  ELEMENT  SIZE 

Figure  2.  Symbolic  texture  Description 
of  raffia. 
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Abstract  representations  of  raffia 
primitives . 


textural  elements  or  primitives  for  a  particular  predominant  element 
size  could  provide  the  average  intensity  and  area  for  that  primitive 
group. 

Implementation  -  In  Fig.  2,  we  can  see  that  al]  of  the 
information  for  a  texture  is  listed  according  to  relative  element 
intensity  and  scan  direction.  For  example,  we  find  that  elements  of 
size  3  and  spacing  8  occur  in  the  vertical  scan  direction;  and  these 
elements  are  dark  in  relation  to  their  vertical  neighbors.  Since  the 
scan  direction  and  relative  element  intensity  of  each  predominant 
element  size  is  known  a  search  for  the  particular  set  of  edge  matches 
exhibiting  these  restrictions  can  be  initiated. 

This  search  proceeds  in  much  the  same  way  that  the  original  edge 
match  search  occurred.  However,  there  are  2  exceptions: 

1)  This  search  is  in  one  scan  direction  for  elements  with  a 
particular  relative  intensity  (dark  or  light)  for  a  restricted 
range  of  distances. 

2)  Instead  of  recording  the  edge  matches  in  an  edge  repetition 
array  the  2  matching  edges  are  marked  in  a  blank  image  in 
locations  corresponding  to  their  locations  in  the  original  edge 
image.  The  interior  points,  points  lying  between  the  2  matching 
edges  along  the  direction  of  scan,  are  also  marked. 

In  Fig.  4  the  marking  scheme  is  illustrated.  Both  edges  are 
marked  by  "E"'s.  The  edge  exhibiting  the  lower  row  location  (and  if 
this  test  fails  the  lower  column  location)  is  marked  in  a  special  way 
as  a  starred  edge.  The  points  between  the  2  edges  along  the  line  of 
scan  are  marked  as  interior  points,  "I". 

The  image  produced  by  this  process  is  the  intermediate  primitive 
image.  This  image  consists  of  zeros  except  where  edge  matches  and 
their  interior  points  are  marked  as  in  Fig.  4. 
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Figure  4 


Figure  5 


00000000 
000  E*  0000 
00010000 
00010000 
00010000 
OOOEOOOO 
00000000 


Non-expanded  primitive  in  vertical 
san  direction. 


00000000 

00000000 

01111100 

00111110 

01111100 

00000000 

00000000 


Expanded  primitive  for  Figure  4. 
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The  next  step  is  the  expansion  of  the  interiors  of  the  primitive 
slices  to  form  a  binary  mask  image  for  the  primitive  set.  Tn  the 
expansion  process  the  original  texture  image  as  well  as  the  edge  and 
direction  images  are  used. 

The  intermediate  primitive  image  is  searched  until  an  "E"  is 
encountered.  Proceeding  in  the  direction  of  scan  the  first  interior 
point,  "I",  is  encountered.  A  "1"  is  placed  in  the  binary  mask  image 
at  this  point  and  the  texture  image  intensity  of  this  point  is  noted. 
This  is  repeated  until  the  edge  mark  at  the  end  of  the  primitive  si  ice 
is  encountered.  The  range  of  intensity  values  for  the  interior  pixels 
of  the  primitive  slice  is  now  known.  At  this  point  the  primitive 
interior  can  be  expanded  to  its  "natural  boundary."  The  natural 
boundary  of  a  primitive  is  made  up  of: 

1)  The  edges  within  the  edge  image; 

2)  pixels  whose  intensities  are  outside  of  the  interior 
intensity  range  by  more  than  10;  and, 

3)  the  boundary  of  the  texture  image  itself. 

Expansion  of  the  primitive  slice  proceeds  outward  from  the 
interior  points  of  the  primitive  slice  in  a  direction  perpendicular  to 
the  direction  of  scan.  Figure  5  shows  the  expanded  primitive  within 
the  primitive  mask  image  corresponding  to  the  primitive  slice  in 
Fig.  4. 

The  expanded  primitives  within  the  mask  image  can  then  be 
analyzed  to  determine  an  average  primitive  intensity,  an  average 
primitive  area  in  pixels  and  from  this  area  the  average  primitive 
dimension  in  the  direction  perpendicular  to  the  line  of  scan. 
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Results 


Results  generated  for  the  raffia  sample  of  Fig.  1  are  given  in 
Fig.  6.  The  dark  primitive  found  in  the  vertical  scan  direction 
corresponds  to  block  B  in  Fig.  3.  The  light  primitive  found  in  the 
vertical  scan  direction  corresponds  to  block  A  in  Fig.  3  and  the  light 
primitive  found  by  the  horizontal  scan  corresponds  to  block  C.  The 
first  primitive  dimension  given  in  each  of  the  primitive  descriptions 
is  the  dimension  along  the  line  of  scan,  while  the  second  dimension  is 
in  a  direction  perpendicular  to  the  scan  line.  Figures  7  through  9 
are  the  primitive  mask  images  corresponding  to  the  primitives  of  type 
A,  B,  and  C  of  Fig.  3,  respectively.  Figure  10  is  composed  of  all  3 
types  of  primitives  with  their  average  intensity  values  kept  constant. 

Work  on  this  program  is  still  in  progress.  It  is  hoped  that 
information  obtained  in  this  way  will  lead  to  the  generation  of  more 
meaningful  texture  descriptions  and  will  be  useful  in  the  texture 
classif ication  process. 
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PRIMITIVE  ANALYSIS  FOR  TEXT.  SUPP12  (THRESH  =10) 

RELATIVE  INTENSITY  IS  DARK  DIRECTION  IS  VERTICAL 
NUMBER  OF  SAMPLES:  108 

AVERAGE  PRIMITIVE  DIMENSIONS  ARE:  (2-00  AND  10. 3*5) 

AVERAGE  PRIMITIVE  SIZE  IN  PIXELS  IS:  (20-30) 

AVERAGE  PRIMITIVE  INTENSITY  IS:  (128.?')) 

PRIMITIVES  REPEAT  AT  ELEMENT  SPACING:  (8-00)  IN  ABOVE  MENTIONED  DIRECTION 

RELATIVE  INTENSITY  IS  LIGHT  DIRECTION  IS  VERTICAL 
NUMBER  OF  SAMPLES:  10B 

AVERAGE  PRIMITIVE  DIMENSIONS  ARE:  (N.QO  AND  ‘I- 33) 

AVERAGE  PRIMITIVE  SIZE  IN  PIXELS  IS:  (3b. <?4) 

AVERAGE  PRIMITIVE  INTENSITY  IS:  (172-35) 

PRIMITIVES  REPEAT  AT  ELEMENT  SPACING:  (8-OD)  IN  ABOVE  MENTIONED  DIRECTION 

RELATIVE  INTENSITY  IS  LIGHT  DIRECTION  IS  HORIZONTAL 
NUMBER  OF  SAMPLES:  b8 

AVERAGE  PRIMITIVE  DIMENSIONS  ARE:  (2-00  AND  7-88) 

AVERAGE  PRIMITIVE  SIZE  IN  PIXELS  IS:  (15.18) 

AVERAGE  PRIMITIVE  INTENSITY  IS:  (I'iC-M?) 

NO  EVIDENCE  OF  PERIODICITY 

Figure  6.  Raffia  primitive  texture 
element  description. 
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Figure  9.  Light  raffia  primitive  found  in 
horizontal  scan. 


Figure  10.  Composite  primitives. 
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2.4  Texture  Edge  Detection 


O.D.  Faugeras  and  H.Y.  Lee 


Background 

The  problem  of  detecting  texture  edges  is  brought  up  in  many 
tasks  of  image  analysis  where  objects  are  present  that  differ  mostly 
in  texture.  Conventional  edge  detection  techniques  based  on  detection 
of  discontinuities  in  intensity  usually  fail  to  locate  boundaries 
between  such  objects.  Thus,  the  need  for  developing  new  techniques 
that  will  perform  better  in  such  cases  comes  out. 

Obviously,  such  techniques  should  be  based  on  some  sort  of 
texture  measurements.  The  quality  that  differentiates  texture 
properties  from  intensity  properties  is  that  the  first  ones  are  local 
that  is  measured  over  some  neighborhood  of  every  pixel  whereas  the 
second  ones  are  pixel  properties. 

This  is  the  first  problem  encountered  when  working  with  texture. 
The  second  problem  comes  in  as  soon  as  we  ask  ourselves  the  auestion 
of  how  to  measure  texture  properties.  To  make  a  long  story  short,  let 
us  say  briefly  that  there  are  two  main  approaches  toward  texture 
analysis.  The  deterministic  approach  considers  textures  as  a  periodic 
or  pseudoperiodic  repetition  of  some  base  pattern  whereas  the 
stochastic  approach  prefers  to  view  texture  fields  as  samples  of 
stochastic  processes.  As  a  consequence  we  can  find  in  the  literature 
two  broad  classes  of  techniques  corresponding  to  these  two  approaches 
which  make  use  of  different  tools:  spectral  analysis,  syntactic 
methods  in  the  deterministic  case,  statistical  methods  in  the 


r 


stochastic  case.  For  a  recent  review  of  the  different  approaches 

toward  texture  analysis  see  fl]. 

Based  on  previous  work  by  Julesz  f 2  ]  and  Pratt  et  a]  .  f31, 
Faugeras  and  Pratt  [4]  proposed  a  technique  of  texture  feature 
extraction  based  on  measurements  of  the  autocorrelation  function  of 
the  texture  field  and  the  first  order  histogram  of  the  decorrelated 
field.  Another  technique  was  proposed  by  Faugeras  (51  based  on  a 
human  visual  model  where  textures  are  analyzed  through  a  bank  of 
nonlinear  channels.  Laws  T6]  expanded  on  both  works  and  proposed  to 
analyze  textures  with  a  set  of  filters  of  small  spatial  extent 
(analogous  to  the  decorrelation  operator  in  f 4 ]  or  the  bandpass 
filters  in  f 5 ] )  and  compute  local  "energy"  values  in  the  output  planes 
after  convolution  (analogous  t  the  L  -norms  in  (5]).  If  N  filters  are 
used  in  the  process,  texture  is  represented  at  every  pixel  by  an 
N-dimensional  vector.  Moreover,  Laws  has  been  able  to  exhibit  a 
limited  set  of  filters  (of  the  order  of  4)  that  seem  to  perform  best 
in  terms  of  classification  accuracy  on  a  fairly  large  set  of  natural 
texture.  Because  of  its  simplicity  and  good  performance,  we  adopted 
Laws'  method  of  texture  feature  extraction  summarized  in  Fig.  1. 

The  problem  of  detecting  texture  edges  is  therefore  equivalent  to 

the  problem  of  detecting  edges  in  the  "energy"  images  I  (k,£)  for 

n 

n=l,...,N.  Edge  detection  methods  fall  in  two  broad  classes. 
Techniques  based  on  enhancement  of  edges  followed  by  thresholding  are 
probably  the  simplest  ones  while  techniques  based  on  local  filtering 
of  an  edge  model  are  slightly  more  complicated.  For  a  good  review  of 
existing  techniques,  see  Davis  (91  and  for  a  quantitative  evaluation 
of  many  methods  falling  in  the  first  class,  see  Abdou  and  Pratt  f7l. 
Two  complications  arise  when  working  with  textures  rather  than 
intensities,  one  of  them  we  already  mentioned  before  it  is  linked  to 
the  fact  that  texture  measures  are  local,  the  other  one  stepping  from 
the  fact  that  we  have  to  deal  with  vector-val ued  functions  rather  than 
scalar-valued  functions  very  much  like  for  the  detection  of  color 
edges.  In  other  words,  we  should  be  able  to  exploit  the 
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I(k,£) 


Figure  1. 


Figure  2 


Texture  feature  extraction  by  parallel  processing  of  an 
image  with  convolution  masks  (k  ,£),...  ,MN  (k,  £)  .  The 
local  energy  measurements  are  given  by 

In(k,£)  =  'y  Jwn  (k-i,  £-j)  Jn  (i,  j)  |  P  n=l ,  .  .  .  ,N 
i  t  j 

where  P  is  a  positive  integer  and  the  window  functions 
w^  may  or  may  not  depend  on  n. 


.  Energy  measures  are  computed  in  neighborhoods  V  and 
V2  of  pixel  (k,£).  The  magnitude  of  their  difference 
is  an  indication  of  the  presence/absence  of  a  texture 
edge  in  the  direction  0  at  that  pixel. 
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cross-correlation  between  the  image  planes  rn^,£)  to  improve  our  edge 
detection  algorithm. 

Texture  Gradient  Computation 

The  method  we  propose  for  detecting  texture  edges  if  patterned 
after  the  methods  for  detecting  intensity  edges  using  directional 
derivatives.  We  can  think  of  this  process  as  acting  either  on  the 
texture  energy  planes  In(k,£)  of  Fig.  1  or  on  the  output  of  the 
convolution  masks  Jn(k,£).  The  idea,  similar  to  what  has  been 
described  in  [8,9]  is  that  if  we  wish  to  detect  a  texture  edge  in  the 
direction  0  (modulo  it)  at  pixel  (k,£),  we  would  want  to  subtract 
energy  measured  computed  in  neighborhoods  and  V2  of  that  pixel  as 
described  in  Fig.  2.  This  is  repeated  for  different  values  of  0,  and 
thus  providing  an  indication  of  the  presence/absence  of  a  texture  edge 
at  every  pixel.  The  actual  detection  is  performed  by  thresholding  or 
local  maximum  selection. 

Now  since  we  have  several  texture  energy  planes  to  work  with,  the 
question  arises  naturally  of  how  to  combine  the  edge  information 
coming  from  those  different  sources.  A  simple  way  of  approaching  this 
problem  would  be  to  select  the  "best"  frames  in  terms  of  some 
criterion  based  on  smoothness  of  edges  and  interimage  edge  element 
correlation.  Another  potential  candidate  is  to  use  some  sort  of 
relaxation  labeling  process  [10,11]  that  would  combine  evidence  from 
the  different  texture  planes  in  a  meaningful  fashion. 

Preliminary  Results 

We  have  applied  some  of  the  ideas  to  the  composite  image  showed 
in  Fig.  3.  It  is  made  of  two  natural  textures,  sand  and  water,  whose 
histograms  have  been  normalized.  This  picture  has  been  processed  by 
the  four  5x5  masks  described  in  Fig.  4.  Vertical  and  horizontal 
texture  edges  have  been  detected  in  all  four  planes  and  results  are 
shown  in  Fig.  5.  Obviously,  masks  #1  and  #2  are  performing  much 
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Figure  6.  Results  obtained  by  combining 
information  obtained  from 
Mask  #1  and  Mask  #2. 
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better  than  the  other  two  masks  and  Fig.  6  shows  the  improvement 
obtained  in  combining  the  results  of  the  first  two  masks. 

As  a  conclusion,  we  believe  that  texture  edges  can  be  accurately 
extracted  by  methods  of  the  type  of  the  one  described  in  this  report. 
We  are  concentrating  now  on  developing  efficient  techniques  of 
combining  the  information  extracted  by  the  different  masks  and  the 
quantitative  analysis  of  the  performances  of  a  texture  edge  detection 
algorithm  based  on  these  ideas. 
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2.5  Application  of  the  General  Linear  Model  to  Binary  Texture 
Synthesis 

D.D.  Garber 


Introduction 


Texture  is  an  important  characteristic  for  the  analysis  of 
images.  Accordingly,  the  study  of  this  image  feature  has  been 
extensive  and  this  has  led  to  the  development  of  models  which  attempt 
to  explain  and  describe  this  quality  of  an  image  field. 
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Earlier  results  [1]  have  indicated  the  usefulness  of  both  Markov 
and  linear  autoregressive  models  in  attempting  to  describe  and 
simulate  natural  textures.  In  the  Markov  model,  a  set  of  generation 
parameters  was  used  to  generate  a  two-dimensional  array  of  binary 
pixels  line  by  line  using  a  set  of  generation  parameters 


3V  <V,,V2 
N+l  L 


•V 


“  P<VN+llVV2' 


>V 


o  > 


where  P(A|B)  represents  the  probability  of  A  given  B.  Each  pixel 

value,  VN+^,  depends  on  the  N  pixels  previous  to  it.  By  allowing 

these  previous  pixels,  V^,  to  be  on  not  only  the  same  image  line  as 

VN+i  but  also  on  lines  above  it  (Figure  1),  horizontal  as  well  as 

vertical  dependence  is  induced  into  the  generating  process.  This 

N+l 

method  of  texture  generation  soon  finds  its  limitations  as  2 

N+l 

generation  parameters  must  be  accounted  for  in  a  binary  case  and  g 
in  a  g  grey-level  case.  It  was  found  that  for  a  fixed  value  of  N, 
better  texture  simulations  could  be  obtained  by  letting  the  pattern  of 
the  IV's  to  become  flexible  and  noncontiguous  (Figure  2). 


This  idea  was  also  used  in  the  linear-model  method  of  texture 
synthesis.  The  linear  model  may  be  expressed  in  normal  linear 
regression  form  as 


Y  =  5t.$+e.  i  =  1,2, ...  ,n  (2) 

i  i  i 

where 


Y  =  V 
i  N+l,i 

i 

< 

M 

> 

'V 

+T 

V2,i 

X2 

Xi  = 

Vi 

X  = 

1 

•  •  •  xz 
_ 1 

1 

5  is  an  (N+l)xl  vector  of  unobservable  parameters  and  e  is  an 

i 

unobservable  random  variable  such  that  Efe.]=0,  The  most  common 
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Figure  1.  Two-dimensional  synthesis 


Figure  2.  One-dimensional  synthesis 


estimator  of  0  is  (X'X)  "''X'Y,  the  least-squares  estimator.  So  given  a 
parent  texture  which  we  desire  to  simulate,  an  estimation  of  the 

A 

texture  model  parameter  0,  0,  may  be  obtained  yielding  the  generation 
model 

Y  =  X0  (3) 

Noise  is  added  to  our  model  during  the  generation  process  which  is 
similar  to  the  noise  we  measure  in  the  original  parent  texture.  It  is 
very  important  to  note  that  the  estimation  of  the  amount  of  noise 
present  in  a  texture  is  usually  obtained  by  measuring  the  amount  of 
variation  in  the  parent  texture  which  is  unexplained  by  our  model . 
Therefore,  as  in  any  modeling  process,  any  shortcomings  or 
inaccuracies  of  the  model  will  appear  to  be  "noise"  (unexplained 
variance)  and  hence  the  amount  of  noise  will  increase.  In  fl)  uniform 
gaussian  noise  was  used  to  drive  the  linear  system  and  the  variance  of 
this  noise  was  estimated  by  measuring  the  variation  of  the  parent 
texture  which  was  unexplained  by  the  linear  model  when  that  linear 
model  was  applied  to  it. 

Choosing  the  V  V s  _to  be  Jji  the  kerne] 

We  will  refer  to  the  '  s  on  which  the  next  pixel,  VN+^ »  depend 

as  the  kernel  of  the  generation  process.  Geometrically  speaking,  the 

\A  's  form  a  kernel  "shape"  which  may  or  may  not  be  spatially 

contiguous.  For  example,  in  Figure  2  a  generating  kernel  shape  is 

shown  where  the  V<-  pixel  depends  on  only  some  of  the  pixels  in  its 

surrounding  neighborhood.  In  this  case,  V  may  be  generated  based  on 

5 

the  values  of  pixels  V-^,  V2 ,  and  but  is  dependent  on  no  other 
pixels  in  the  neighborhood.  This  does  not  imply  that  is  not 
related  or  correlated  with  its  other  neighbors.  In  fact,  the 
relationships  between  and  V  will  determine  other 

interrelationships . 
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A  non-contiguous  neighborhood  of  V^'s  is  used  as  it  allows  a  more 
parsimonious  model  for  texture  generation  to  be  chosen.  An  analogy  in 
the  simple  linear  regression  case  would  be  the  dropping  of  independent 
variables  which  do  not  contribute  to  the  prediction  or  estimation  of 
the  dependent  variable.  In  texture  generation  this  allows  the  model 
to  be  estimated  by  fewer  parameters  and  makes  the  generation-synthesis 
process  more  efficient.  It  also  relieves  us  from  the  complex 
numerical  problems  of  analyzing  and  inverting  matrices  of  unwieldy 
size,  a  necessary  step  in  linear  model  parameter  estimation.  We 

would,  for  example,  not  expect  our  VN+^  pixel  to  depend  on  a  pixel 
where  the  spatial  separation  between  VN+1  and  is  large.  If  that 
distance  is  small,  however,  we  would  expect  a  large  dependence. 

Methods  for  choosing  the  proper  independent  variables  (V^'s)  to 
be  included  in  the  generation  process  have  been  developed  over  many 
years  [3,4,5]  for  other  applications.  If  we  limit  ourselves  to  a 
finite  neighborhood  of  M  variables  we  would  wish  to  choose  the  best 
subset  of  N  variables,  where  N<M  and  is  more  manageable  in  size. 
Evaluating  such  subsets  and  their  corresponding  models  requires  a 
criterion  function.  The  most  commonly  used  is  mean-square  error.  The 
problem  then  is  to  choose  a  subset  of  V^'s  of  size  N  from  a  set  of 
V.'s  of  size  M  (M>N)  such  that  the  model  Y=Xg+e  employing  those  V^'s 
produces  the  minimum  mean-square  error  when  compared  to  all  other 
possible  models  containing  N  V.'s.  The  only  possible  way  to  do  this 
is  to  evaluate  all  possible  {^)  combinations  of  points.  Even  for 
small  values  of  M  and  N  this  would  require  enormous  computation  power. 
For  example,  for  M=40  and  N=12,  over  5.5  billion  models  would  have  to 
be  evaluated!  Therefore  a  sub-optimal  approach  which  yields  a  good 
solution  but  not  necessarily  the  best  solution,  must  be  used. 

Correlation  and  Partial  Correlation 

We  may  define  the  mean  and  variance  of  a  variable  X  to  be 
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Ux  =  E[x] 


M) 


and 


a l  =  E[ (X-y)2] 


(5) 


Similarly  we  may  define  the  covariance  of  two  variables  X^  and  X^  as 


Yy  y 
X1X2 


=  E[ (X1-yx  ) (X2-Vx  > 1 
EC (X1X2)-yx^yx^] 


(6) 


And  their  correlation  coefficient  as 


X,  X 


12 


'X1X2 

°Xl\ 


H) 


Using  the  above  definitions  we  may  then  define  the  partial  correlation 
coefficient  of  variables  X ^  and  X2  after  both  have  been  adjusted  for 
X3  as 


P12  •  3 


p12  tp13*  (p23) 


( R ) 


Each  of  the  above  parameters  has  a  corresponding  statistic  (estimate 
of  a  parameter  of  a  population  given  an  observable  sample  of  that, 
population)  which  is  chosen  to  meet  some  desirable  set  of  a  criteria 
such  as  sufficiency,  consistency  and  unbias  of  estimate  under  certain 
conditions.  These  would  be 
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Higher  order  partial  correlations  maybe  found  by  extension  of  the 
above  [7 ] . 

The  Forward  Selection  P rocedu re 

One  approach  to  finding  a  good  subset  of  independent  variables  is 

known  as  the  forward  selection  procedure.  This  method  inserts  N 

variables  into  the  model  one-at-a-time .  The  order  of  insertion  is 

determined  by  using  the  partial-correlation  coefficient  as  a  measure 

of  the  importance  of  variables  not  yet  in  the  equation.  The  basic 

procedure  is  as  follows.  First  we  select  the  V.  most  correlated  with 

^  1 

VN+i  and  we  find  the  first-order  linear  regression  equation 

V  =aV-+b.  We  next  find  the  partial  correlation  coefficient  of  V 
N+l  H  3 

(for  all  j=i)  and  VN+1  (after  allowance  for  V^) .  Mathematically  this 
is  equivalent  to  finding  the  correlation  between  the  residuals  from 

A 

the  regression  vN+i=aiVi  +a2  an<3  residuals  from  another  regression 

Vi  +b2  (which  we  liave  not  actually  performed).  The  with  the 

highest  partial  correlation  coefficient  with  VM, . ,V.  is  now  selected 

/\  £M"r  ±  1 2 

and  a  second  equation  V.,,  ,=c.V.  +c_V.  +c_  is  fitted  to  the  data.  This 

N+l  1  i.  2  i2  3 

process  continues.  After  V.  ,V.  ,...,V.  are  in  the  regression  the 

■*-1  12  xq 

partial  correlation  coefficients  are  the  correlations  between  the 

A 

residuals  from  the  regression  V  ,=f(V.  ,V.  ,...,Vi  )  and  the 

A  INi-L 

residuals  from  a  regression  V^.=f.(Vi  .V^  ,...,V^  )  (j>q).  The  Vi . 

«•  j  3  1  2  q  3 

with  the  highest  partial  correlation  coefficient  is  now  selected  for 

entry  into  the  linear  model.  The  process  is  continued  until  N  ' s 

are  entered  into  the  model. 

The  final  N  variables  chosen  by  a  forward  selection  procedure  is 
not  guaranteed  to  be  the  optimal  set  but  given  the  logistics  of  the 
selection  procedure,  the  solution  obtained  is  usually  close  to 
oct imal . 


Numerical  Implementation 


One  of  the  most  common  procedures  for  implementing  the  forward 
selection  process  numerically  utilizes  Doolittle  decomposition  T5]. 
The  Doolittle  decomposition  may  be  used  to  find  the  inverse  of  the 
correlation  matrix,  R,  and  the  estimates  of  linear  model  coefficients 
as  each  variable  is  entered  in  the  model.  The  correlation  matrix  is 
factored  into  the  product  of  two  triangular  matrices 


Rp  =  LpUp  (14) 

where  Lp  is  lower  triangular  and  Up  is  upper  triangular  with  ones  on 
the  diagonal.  Partial  autocor relat ion  coefficients  may  be  obtained 
easily  from  elements  of  matrix  during  its  decomposition  at  each  step. 
For  further  details  and  examples  see  Beyer  M]  and  Draper  and  Smith 
[3]. 

Choosing  N 

In  the  practical  case  N  is  often  chosen  by  computational  limits 

imposed  by  finite  processing  capability  or  finite  memory.  The  ideal 

of  parsimony  would  also  urge  the  selection  of  the  smallest  N  possible. 

In  the  most  general  texture  synthesis  (Markov)  model  which  utilizes 

N 

N-grams  [1]  we  find  that  as  many  as  g  storage  locations  are  required 

in  order  to  collect  the  data  needed  to  synthesize  a  texture  from  a 

parent  where  g  is  the  number  of  grey  levels  in  the  texture.  This 

approach  calls  for  N  to  be  less  than  17  for  g=2  (a  binary  image)  if  we 
17 

have  only  2  =131,072  stora'  locations  in  memory.  Approaches  to 

"stretch"  this  limitation  have  been  investigated  and  will  be  discussed 
in  later  papers.  If  we  have  8  grey  levels  then  N  must  be  less  than  or 

5 

equal  to  5  as  8  =32,768.  Thus,  in  synthesizing  textures  using  an 
N-gram  approach,  processor  storage  capability  is  the  major  limiting 
factor  . 
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The  general  linear  model  approach  to  synthesizing  textures  was 
developed  to  overcome  this  limitation  and  produce  a  simpler  model. 
However,  as  was  discussed  earlier,  model  parameter  estimation  requires 
the  inversion  of  an  (N+l)  matrix  where  N  is  the  number  of  variables 
in  the  model  (1 ) . 

Application  to  Binary  Textures 

The  linear  autoregressive  approach  may  be  used  to  generate  binary 
textures  as  well  as  continuous-tone  textures.  Applying  such  a  model 
to  binary  textures  corresponds  closely  to  using  a  linear  discriminant 
function  in  a  two-class  discrimination  problem.  A  decision  (white  or 
black)  is  made  based  on  a  linear  combination  of  previous  measurements. 
Two  parent  textures  water  (Figure  3a)  and  raffia  (Figure  4a!  were 
chosen  and  estimates  of  parameters  of  the  linear  models  were 
determined.  Figure  3b  (N=56)  and  4b  (N=70)  show  the  binary  texture 
generated  using  the  resulting  linear  model.  The  noise  used  in  the 
generation  process  was  normally  distributed,  although  better  choices 
in  this  case  could  possibly  be  found.  At  each  pixel  generation,  the 
value  of  the  next  pixel  (obtained  by  taking  a  linear  combination  of 
previous  binary  pixels  in  its  neighborhood  (Figure  2)  plus  noise!  was 
quantized  ( thresholded )  to  a  binary  value  (white  or  black!  before 
generating  the  next  pixel.  Texture  synthesis  using  the  more  complex 
Markov  model  are  also  shown  in  figures  3c  and  4c.  Additional  methods 
of  synthesizing  binary  textures  using  the  ideas  of  Cox  (8]  are 
currently  being  investigated. 

Application  to  Discrimination  and  Identi f ication  Textures 

By  making  assumptions  concerning  the  distribution  of  the  error 
term,  T,  in  the  linear  model  Y=£X+t  which  is  fit  to  the  data, 
probability  distributions  nay  be  derived  for  the  statistics  of  the 
regression  (autoregressive)  model.  It  has  been  shown  by  Graybill  r2i 
and  Box  and  Jenkins  [5]  that  the  statistical  estimate  of  i,  t  has  mean 
5  and  covariance 


77 


'  Jjaaix 


8 


(15) 


A 

cov [8] 


-1  2r 
n  a£r 


-1 


where  r  is  the  auto  covariance  matrix  of  the  N  successive  values  of 
the  time  series  process.  To  derive  (15)  it  is  necessary  to  assume 
that  the  elements  of  e  are  independent  with  constant  variance  a£. 
Estimates  of  the  variances  and  covariances  are  obtained  by 
substituting  estimates  for  the  parameters  in  (15)  yielding 


A  ^  ^  A  a  A  « 

covie]  =  n~1azT~1 
e 

~  —It 
r  =  n  •L(Xi  X) 

where 

Z-N+l 
Z-N+2 

X  = 

.  Zn-N 


Z-N+2  Z-N+3  * ‘ *  Z0 

Z~N+3  Z-N+4  *  *  *  Z1 

* 

• 

•  •  •  z  , 

n-1 


(lfi) 


(17) 


and 


A 

a 


2 


n  _  -I  n  « 

{  E  [z+-n  x  E  z  1  } 
t=l  T  u=  1  u 


For  a  two-dimensional  time  series  such  as  a  texture,  the  z^  are 
obtained  by  passing  the  chosen  generation  kernel  containing  N  points 
over  the  texture  of  interest  n  and  measuring  the  values,  V.  ,  at  each 

i 

of  the  the  n  positions. 
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Given  this  fact,  it  is  easy  to  test  various  hypotheses  about  g 
and  confidence  limits  may  be  obtained  on  its  individual  elements  if  a 
distribution  of  the  errors  e  ^  is  assumed.  Sub-model  tests 
incorporating  subsets  of  the  g vector  may  also  be  tested.  Basically, 
at  each  step  of  the  forward  selection  procedure,  a  hypothesis  testing 

A 

whether  6.  =0  is  performed  where  g .  is  the  estimated  linear  coefficient 
corresponding  to  the  newly-entered  variable  V  at  that  step.  This  is 
accomplished  by  observing  the  test  statistic 


incremental  sum  of  squares  due  to  V./l 
F  =  _ i _  (IP) 

sum  of  squares  residual/ (N-i-1) 

This  statistics  has  an  F  distribution  with  degrees  of  freedom  1  and 
(N-i-1)  when  the  errors,  e^,  of  the  model  are  normally  and 
independently  distributed.  When  all  V^'s  not  included  in  the  model  at 
step  i  fail  to  meet  the  test  of  significance  the  process  of  adding 
variables  to  the  model  may  be  stopped.  Unfortunately,  halting  the 
process  at  this  point  does  not  insure  statistical  of  all  variables 
entered  to  that  point  nor  does  it  insure  non-significance  of  variables 
not  entered.  This  problem,  as  that  of  finding  the  optimal  set  of  N 
variables,  may  be  solved  only  by  considering  all  possible  models. 

The  distribution  of  the  3^'s  can  be  used  to  test  the  equivalence 
of  two  linear  models  and  therefore,  the  equivalence  of  the  textures 
which  those  models  attempt  to  explain.  Application  of  this  concept  to 
identification  and  discrimination  of  texture  fields  will  be  done  in 
future  work. 

Conclusion 


The  application  of  the  concept  described  in  this  paper  led  to  a 
good  simulation  of  a  variety  of  textures  using  both  the  Markov  and 
linear  autoregressive  model  T 1 ] .  An  approach  for  applying  the  linear 
model  to  the  problem  of  texture  identification  and  discrimination  was 
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proposed  and  is  currently  being  investigated. 
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2.6  Algebraic  Reconstruction  Techniques  for  Texture  Synthesis 
O.D.  Faugeras  and  D.D.  Garber 


Introduction 

Julesz's  conjecture  [1]  that  second-order  statistics  are 
sufficient  in  terms  of  human  visual  texture  discrimination  provides  a 
useful  bound  on  the  amount  of  information  necessary  to  reconstruct  a 
texture  field.  Although  counterexamples  to  that  conjecture  have 
recently  been  found  [2,3],  it  still  seems  to  be  a  fairly  accurate 
first-order  approximation.  Examples  of  the  use  of  that  upper  bound 
for  texture  analysis  can  be  found  in  [4,5].  Therefore  it  is  very 
tempting  to  use  it  for  synthesizing  natural  texture  fields.  Previous 
approaches  to  texture  synthesis  have  used  the  nth-order  Markov  model 
approach  [3,6,7],  the  autoregressive  model  approach  [6,8,91  or  the 
syntactic  approach  [10]. 

The  purpose  of  this  paper  is  to  show  that  if  we  use  the 
generation  coefficient  approach  (to  be  defined  below)  to  synthesize 
textures,  even  though  we  limit  our  knowledge  of  the  original  random 
field  to  second-order  statistics  (and  in  fact  to  kth-order  statistics, 
k  2  2)  we  are  confronted  to  the  task  of  "inventing"  higher-order 
statistics.  We  will  demonstrate  how  this  can  be  done  using  Algebraic 
Reconstruction  Techniques  (ART) . 

Problem  Definition 


The  stochastic  approach  toward  Texture  Analysis  considers  texture 
fields  as  samples  of  2-D  stochastic  fields.  Assuming  that  we  are 
dealing  with  sampled  and  quantized  imagery,  let  f(n^,n2)  denote  the 
random  field.  n^  and  n 2  are  integers  representing  coordinates  of 
points  in  the  plane.  Let  n  be  the  vector  of  coordinates  n^  and  n^» 
Second-order  statistics  are  given  by  the  set  of  second-order  joint 


* 

L  — 
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density  functions 


n 


,S(VX2> 


(1) 


for  all  possible  vectors  n  and  in  ,  where  x^  and  x^  are  the  values  of 
the  random  variables  f(n)  and  f(in),  respectively.  We  will  also  assume 
for  simplicity  that  the  random  field  is  homogeneous,  that  is, 
invariant  through  translations: 

p->  ->  =  p->-  -+■  (?) 

n ,  m  rn-m 

Suppose  that  we  are  given  the  set  of  joint  density  functions  p+  for 
all  vectors  r  and  that  we  want  to  synthesize  a  texture  field  with  the 
same  second-order  statistics.  Since  we  have  to  work  with  finite 
images,  we  will  assume  that  1  _<  <  N. 

There  are  many  ways  of  carrying  out  such  a  synthesis.  If  we 
stick  to  a  television  type  of  scanning  (left  to  right,  top  to  bottom), 
we  will  first  generate  f(l,l)  using  first-order  statistics,  then 
f  ( 1 , 2 )  using  second-order  statistics,  f(l,3)  using  third-order,  etc... 
Even  if  we  assume  finite  memory,  namely  that  intensity  at  joint  (i,j) 
depends  only  on  intensities  at  points  located  in  some  finite 
neighborhood,  we  see  that  second-order  statistics  are  not  sufficient 
to  generate  the  process.  Indeed,  assuming  a  memory  of  order  P,  the 
intensity  at  joint  (i,j)  is  computed  using  the  conditional  probability 
or  generation  coefficient 


p  (x ,  . 
^  il 


1J1 


r  X  .  .  ) 

Iplp 


(3) 


which  involves  (P+1 ) th-order  joint  density  functions.  Therefore  if  P 
is  larger  than  1,  we  have  to  "invent"  those  higher-order  densities. 


i 
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Another  way  of  implementing  the  generating  process  would  be  to 

draw,  at  random,  pairs  of  coordinates  (i,j)  and,  again  assuming  finite 

memory,  to  examine  the  Q  points  (0  Q  <  P)  already  generated  within 

the  neighborhood  of  point  (i,j).  Intensity  at  that  point  wou]d  then 

be  drawn  conditionally  to  the  Q  already-drawn  intensities.  But  again 

we  would  have  to  "invent"  higher-order  joint  densities. 

Mathematically  the  problem  can  then  be  stated  as  follows:  given  M 

random  variables  X, ,...,Xw  such  that  for  every  pair  (X.,X.),  1  <  i,j  < 

1  M  1  j  —  — 

M  and  i  /  j ,  we  know  the  joint  density  function  p(x^,x.),  find  a 

function  p(x, ,...,xw)  which  satisfies 
1  M 


p(xlf . . . ,xM)  >  0 


for  all  x^, . . . ,xM 


xR,  K^i  and  K/j 


P(xlf... ,xM)  =  P ( x^ , x j ) 


The  p(x^,Xj)'s  are  sometimes  called  the  marginals  of  p(x^, . . . ,xM) . 
Assuming  quantization  with  K  levels,  the  K14  unknowns  p(x1,...,xM)  can 
be  stacked  up  in  a  vector  ]5  and  conditions  (4)  and  (5)  correspond  to  ^ 
linear  programming  problem: 

/  Ap  =  in  ( f  ) 


p  >  0 


where  vector  in  is  obtained  from  the  functions  p(x.  ,x.)  and  matrix  A  of 

w  9  M  ^  ^ 

size  ((”)K  )xK  contains  only  ones  and  zeroes.  For  any  reasonable 
values  of  M  and  K,  this  is  a  set  of  linear  equations  and  inequalities 
of  fairly  large  dimension  and  usual  resolution  techniques  like  the 
simplex  [3]  become  very  quickly  limited. 


Solution  through  Algebraic  Reconstruction  Techniques  (ART) 

ART  was  introduced  by  Gordon,  Bender  and  Herman  I'll]  for  solving 
the  problem  of  three-dimensional  reconstruction  from  projections  in 
electron  microscopy  and  radiology.  This  is  a  deconvolution  problem 
where  an  estimate  of  a  function  in  a  higher-dimensional  space  is 
deconvolved  from  its  experimentally  measured  projections  to  a 
lower-dimensional  space.  For  an  excellent  review  of  those  techniques 
see  Gordon  [12]. 


The  problem  stated  in  equations  (4)  and  (5)  or  (6)  and  ('>)  is 
precisely  a  problem  where  the  projections  are  the  second-order  joint 
density  functions.  ART  is  therefore  directly  applicable.  The  simple 
intuitive  basis  is  that  each  projected  density  is  thrown  back  across 
the  higher-dimensional  region  from  whence  it  came,  with  repeated 
corrections  to  bring  each  projection  of  the  estimate  info  agreement 
with  the  corresponding  measured  projection. 

Formally,  we  use  an  iterative  scheme  defined  by 


^(q+D 


= 


(q) 


P  J  /  •  •  •  /  P  (X^****/  +tc  (^i  f  •  •  ft  r  x^) 


for  all  values  of  x.  ,.  .  .  ,x.. 

1  M 
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for  q  =  0,1,... 


the  correction  term  c  ^*(x^,...,x  )  is  given  by: 


(q> 


c (Xj_,  .  .  .  ,xM)  2  ^  ^  (p(xifXj)-p(q)  (x^Xj)  ) 

'2'  i=l  j=i+l 


where  p { x 
original 


i 


»x.)  is 

texture 


the  actual 
field  and 


marginal  measured, 

p(q^(x,,x,l  is  the 
1  3' 


for  example,  from  an 
marginal  corresponding 
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to  the  reconstructed  density  at  iteration  q. 


We  may  express  this  in  words  as  follows.  The  iterative  process 

may  be  started  with  all  reconstruction  elements  set  to  a  constant 

(p*°* (xlf . . . ,x  )  =  for  all (x1# . . . ,xM) .  In  each  iteration  the  sum  of 

the  dif ferencesKbetween  the  actual  and  the  reconstructed  marginals  is 

M—  2 

computed  and  evenly  divided  amongst  the  K  reconstruction  elements. 
If  the  correction  is  negative,  it  may  happen  that  the  calculated 
density  becomes  negative  at  some  points.  This  problem  can  be 
alleviated  by  using  a  modified  iteration  scheme  defined  by 


(q+1) 


(x, 


'V 


=  Max{  0  ,p  ^  (x 


,xM)+tC (q) (X1' 


’XM)} 


(10) 


therefore  guaranteeing  p^+1  0  (constrained  ART  [12]).  It  is  of 

course  necessary  to  determine  when  an  iterative  algorithm  has 
converged  to  a  solution  which  is  optimal  according  to  some  criterion. 
This  is  in  turn  related  to  the  problem  of  finding  the  optimal  value  t 
of  t.  Various  criteria  for  convergence  have  been  devised  [12].  For 
simplicity,  we  chose  the  mean-square  error  between  the  measured  and 
calculated  marginals: 


sq  -  I|5(q)!l2  =  llS-S<q)ll2  <”> 

where  (|.||  is  the  usual  euclidean  norm  and  m  is  defined  in  equation 
(6).  The  derivation  is  given  in  the  Appendix. 


A  dual  approach  that  we  explored  also  is  based  upon  the  analysis 
of  equation  (5)  in  the  Fourier  domain.  It  can  easily  be  shown  that 
the  initial  problem  stated  by  equations  (4)  and  (5)  is  equivalent  to 
an  interpolation  problem  in  the  Fourier  domain.  The  major  drawback  of 
this  approach  is  the  difficulty  to  ensure  the  positivity  of  the 
inverse  Fourier  transform  of  the  interpolated  function.  Therefore  we 
did  not  pursue  in  that  direction  even  though  it  may  be  the  case  that 
"good"  interpolating  functions  will  alleviate  that  problem. 
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The  basic  philosophy  of  the  two  approaches  we  just  discussed  is 

that  Mth-order  joint  density  functions  are  "invented”  to  satisfy 

exactly  the  constraints  stated  in  equation  (5).  Their  obvious 

disadvantage  is  the  high  dimensionality  (KM  for  an  Mth-order  joint 

density  function)  of  the  data  that  is  to  be  stored  compared  with  the 

.  M  2 

usually  lower  dimensionality  ((^)K  for  the  second-order  joint  density 
functions)  of  the  data  that  is  effectively  used.  An  obvious  way  of 
minimizing  this  problem  is  to  relax  the  constraints  of  equation  (5) 
and  be  content  with  functions  p(x1#...,x  )  which  satisfy  on] y  a  subset 
of  these  constraints.  Chow's  expansion  ([13],  Chapter  4)  paves  the 
way  to  such  an  approximation  with  the  advantage,  when  compared  with 
other  expansions  such  as  the  Rademacher-Walsh  and  the 
Bahadu-Lazarsfeld ,  that  it  always  provides  a  full-fledged  joint 
density  function,  that  is  positive  and  summing  to  ].  Chow's  expansion 
can  be  written  as 

p<xl . V  =  P(x1)p(x2|x1)...p(xH|x1,x2 . XM_1)  ,1?' 

This  is  an  exact  representation  and  as  proposed  in  r 1 4  ]  ,  product 
approximations  can  be  found  that  use  only  second-order  densities. 
Their  general  definition  is 

M 

p(xl . V  =  p<xm.  IV  ... 

i=l  i  k(i) 

where  (m  ^,  .  .  .  ,m  ^  is  some  permutation  of  (1,...,M),  m  =0  and  pfx^Jx  ) 
=  p(x.  ).  It  can  be  easily  shown  that  p  is  a  valid  density  function 
and  that  if  k(i)  equals  0  for  only  one  value  of  i,  equation  (5)  holds 
for  M-l  pairs  of  indexes  i  and  j.  Unfortunately  we  have  not  been  able 
to  obtain  good  synthesis  using  this  approach  at  the  present  time. 


0  <  j  (i)  <  i  H31 
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Results  and  Conclusions 

The  iterative  process  defined  in  (8)  may  be  halted  at  any  number 
of  iterations,  q,  and  a  texture  may  be  generated  using  the  value  of 
at  that  point.  However,  it  should  be  kept  in  mind  that  the  success  of 
a  texture  synthesis  depends  on  making  the  error  c^  as  small  as 
possible  and  the  texture  generation  process  is  sensitive  to  this 
error.  It  has  also  been  found  by  experimentation  that  the  p  contains 
many  values  which  are  set  to  zero  by  implementation  of  constrained  ART 
(10).  This  tends  to  cause  the  Markov  texture  generating  process  to 
become  absorbing,  which  causes  patches  of  white  and  black  or  streaks 
and  lines  to  be  generated.  This  is  eliminated  by  setting  those  values 
which  are  zero  to  a  small  non-zero  value,  e  ,  in  the  generation 
process . 

Using  the  above  concepts,  texture  simulations  of  the  binary 
textures  water  (Figure  la)  and  raffia  (Figure  lb)  were  generated 
(Figures  2a, b).  Synthesized  textures  employing  actual  kth-order 
statistics  (K=14)  were  also  generated  (Figures  3a, b).  Simulations  of 
other  textures  employing  kth-order  statistics  and  details  of  the 
texture  generation  method  used  are  contained  in  rs]. 

Append i x 

Derivation  of  t 

Rewriting  equation  (8)  in  vector  form,  we  get: 

^<q+i>  =  ~(q)  +  ta(q)  (An 

Multiplying  both  sides  of  equation  ( A  .1  )  with  matrix  A  (equation  (M  ) 
we  obtain 


~  (q+D 

m 


^(q> 

=  m  ^ 


ta(q) 


(A?> 
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where 


«  Ae(q)  (A3) 

subtracting  the  actual  marginal  vector  m  from  both  sides  of  equation 
(A2)  yields 

||e(q+1)||2  =  ||e(q)-t3(q)||2  =  t2||3(q)||2-2t2(q)-e(q)+||e(q,||2  (A4) 

Therefore  the  error  eq+i  at  iteration  q+1  is  minimized  for 


.^(q) 
tq=  l|3(q>||| 


(A5) 


where  *  denotes  the  inner  product. 
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2.7  Autoreg res i ve  Modeling  with  Conditional  Expectations  for  Texture 
Synthesis 

O.D.  Faugeras 


Introduction 


Julesz'  conjecture  that  second  order  statistics  are  sufficient  in 
terms  of  human  visual  texture  discrimination  ill  provides  a  useful 
bound  on  the  amount  of  information  needed  to  reconstruct  a  texture 
viewed  as  a  realization  of  a  random  process.  Even  though  this 
conjecture  has  recently  been  disproved  [2,3],  it  nonetheless  provides 
an  excellent  first  order  description  of  human  texture  perception.  For 
an  investigation  of  the  tightness  of  that  bound  see  f4]  and  for  an 
application  to  texture  analysis  see  [5],  We  show  in  a  companion  paper 
[6]  that  even  with  this  upper  bound  an  approach  to  synthesis  based 
upon  the  use  of  high  order  density  functions  is  rapidly  limited  by  the 
size  of  the  data  to  be  handled.  On  the  other  hand,  an  approach  based 
only  upon  the  linear  prediction  of  the  intensity  value  at  one  pixel 
from  neighboring  pixels  intensities,  although  quite  simple,  takes  only 
into  account  the  autocorrelation  function  of  the  texture.  Early  work 
on  the  subject  by  McCormick  and  Jayaramurthy  [7]  made  use  of  a 
one-dimensional  autoregressive  model.  Deguchi  *nd  Morishita  C81,  Tou 
et  al.  [9]  and  Tau  and  Chang  [10]  also  used  a  similar  technique. 
Garber  [11]  extended  these  ideas  to  a  two-dimensional  model. 
Obviously,  except  for  binary  fields  or  gaussian  fields  this  is  a  lot 
less  information  than  what  is  contained  in  second  order  statistics. 
These  two  reasons  prompted  us  to  explore  an  intermediate  approach  in 
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which  we  linearly  predict  the  intensity  value  of  one  pixel  from  the  , 
expected  values  of  that  intensity  given  the  intensity  values  of  sets  j 
of  neighboring  pixels.  The  advantages  of  this  approach  are  that  high 
order  statistics  can  be  taken  into  account  in  the  synthesis  process  in 
a  controlled  manner  while  keeping  the  simplicity  of  the  linear 
prediction  techniques. 


We  will  describe  in  this  paper  a  method  suitable  for  analyzing 
and  synthesizing  so  called  "natural"  texture  fields  where  the  word 
natural  means  that  there  is  a  sufficient  amount  of  randomness  so  that 
a  modelization  based  upon  stochastic  fields  theory  is  not  completely 
preposterous.  We  therefore  eliminate  from  our  scope  purely  repetitive 
patterns  probably  best  described  by  a  structural  approach. 


The  Model 


Let  )^(n  be  a  discrete  stochastic  field  where  integers  m  and  n 
are  the  space  coordinates  and  xm  n  is  the  intensity  at  point  (m,n) . 
For  convenience  in  the  notations  we  will  use  only  one  index  k  to 
denote  the  pair  (m,n) ,  k  may  therefore  be  thought  of  as  a  vector  of 
coordinates  (m,n) . 

Linear  prediction  techniques  attempt  to  approximate  the  random 
variable  xk  by  a  weighted  sum  of  the  random  variables  at  neighboring 
points. 


°iXi 


(1) 


The  a^'s  are  then  adjusted  to  minimize  the  variance  of  the  error 

Et<xk-*k>2} 

where  £{•}  denotes  the  expected  value.  This  is  achieved  simply  by 
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solving  a  set  of  linear  equations 

Aa  =  t 

->•  T 

where  matrix  A  is  NxN  symmetric  with  a^  -  Efx^Xj},  a  =  fai»...raNl 
and  b  »  fE  (x^ ,  x  , . . .  ,E{xk,xN}]  T.  Therefore,  the  o^'s  contain  only 
information  about  the  autocorrelation  function  of  the  stochastic 
field. 


Our  attempt  to  use  more  than  second  order  moments  while  keeping 

the  simplicity  of  the  linear  prediction  approach  is  based  on  using  the 

expectations  of  the  value  at  point  k  conditioned  on  the  values  at 

neighboring  points.  Very  simply,  if  we  know  the  random  variables 

x.  ,...,x.  ,  the  best  prediction  of  x.  is 
X1  Xq  K 

E{xklxi1,*‘*/Xi  } 

1  q 

In  other  words  E{x.|x.  X|  }  is  the  function  of  Xj,  ,...,x^  which 

*  i  p  p 

provides  the  best  approximation  of  x^  in  the  mean-square  sense  T 1 2 1 . 
Thus,  we  define 

K 

“iE{xklyi) 


as  as  estimate  of  x,,  where  y.  denotes 

k  1 1 


the  q-tuple 


of 


(xi1"*-xiq) 

random  variables  at  nearby  points.  Notice  that  computing  Eix^Jy^} 
implies  that  we  know  the  q+lst  order  joint  density  function  pix^^). 
The  choice  of  q  and  K  is  governed  by  our  willingness  to  achieve  a  more 


accurate  prediction.  Once  this  choice  has  been  made, 
computed  so  that  they  minimize  the  error  variance 


the 


ai '  s 


are 


E{  (x^-x^)2} 


The  mathematics  of  the  derivation  are  essentially  the  same  as  in  the 


simple  case  of  the  linear  prediction  technique.  The  details  are  given 
in  Appendix  A.  We  again  have  to  solve  a  set  of  linear  equations. 


Aa  =  S 


(4) 


A  is  a  KxK  symmetric  matrix  whose  current  element  a. .  is  given  by 


lij  ~  E{E^kIyi>E{xk|y.}} 

=  S-  E{xklyi}E(xkiyj)p(yi'yj) 


(5) 


1  D 


vector  b  's  ith  coordinate  is 


b±  =  EtxkEtxklyi>} 


£ 

xk'yi 


xkE{xklyi}p(xk,Yi) 


(«) 


p(yi»yj)  is  the  (2q)th  order  joint  density  function  of 


the  random 


variables  (xj  )  and  (x-i,  , . .  . , x.:  )  while  p(xwy.)  is  the  q+lst 

xq  Jq  k.  i 

order  joint  density  function  of  the  random  variables  x,  and 


(X^  f  ...  X  4  ). 

1  q 


Although  everything  we  said  is  valid  in  the  continuous  case  where 
the  random  variables  xk  can  take  any  values  we  deal  in  practice  with 
quantized  values.  We  will  therefore  assume  in  what  follows  that  the 
xk's  take  only  G  integer  values  (grey  levels)  which  we  can  assume 
without  restriction  of  generality  to  be  Ofl,...,G-l. 


Complexity,  Synthesis  Method 

The  model  contained  in  Eq.  (3)  is  made  of  two  different  parts. 
First  the  conditional  expectations  E{xk|y,}  have  to  be  computed  from 
some  real  data  by  estimating  the  q+lst  order  joint  density  function 
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p<VV 


VJ 

Ej^lyj)  -  2  xkp(xk|yi) 


*  x  =0  ~ 
k 

There  are  values  to  compute  and  store.  Second,  numbers  are 
computed  by  solving  Eq.  (4).  This  implies  estimating  the  (2q)th  order 
joint  density  functions.  p(y^»yj)  for  i  =  1,...,K;  j  =  i+l,...,K  each 
one  being  represented  by  G^  values.  The  texture  is  then  represented 
by  the  KxG^  numbers  E{x,  |y.}  and  the  K  numbers  a.. 

We  adopted  for  simplicity  a  synthesis  method  based  on  a  top  to 
bottom,  left  to  right  scanning  of  the  image.  Therefore,  the  actual 
neighborhood  used  to  predict  the  value  x^  at  point  k  is  as  indicated 
in  Fig.  1.  At  this  point  it  may  be  worthwhile  to  discuss  problems 
which  we  encountered  when  we  tried  to  implement  these  ideas.  The 
first  one  is  related  to  the  model  that  we  used.  In  order  to  generate 
an  actual  texture  field  (for  example  a  512x512  image)  the  process 
described  by  Eq.  (3)  has  to  be  seeded.  For  points  in  the  image 
located  close  to  the  first  row  or  the  first  and  last  columns  intensity 
values  have  to  be  set  at  neighboring  pixels  lying  outside  the  image. 
These  intensities  are  usually  taken  as  values  obtained  from  a 
pseudo-random  number  generator.  But  even  with  such  random  initial 
conditions,  the  process  described  by  Eq.  (3)  converges  quickly  toward 
a  deterministic  steady  state  solution.  To  avoid  that  problem  (also 
encountered  in  [5])  we  can  rewrite  Eq .  (3)  as 


A 

*k  =  S  aiE{xklyi}  +  PUk 


where  the  process  is  a  white  process  composed  of  zero  mean,  unit 
variance  equid istr ibuted  random  variables  independent  of  the  random 
field  x  .  Because  of  these  hypothesis  nothing  changes  in  the 


97 


o 

o 

o 

o 

o 

o 

o 

o 

o 

• 

• 

• 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

X 

— 

-> 

*k 

>  direction  of  scan 


Figure  1.  Neighborhood  used  to  predict  the  value 


\ 


determination  of  the  numbers  ou  but  the  parameter  can  be  adjusted  so 

that  the  variance  of  x,  is  the  same  as  the  variance  of  x,  (see 

k  k 

Appendix  B)  . 


The  second  problem  that  we  encountered  is  related  to  an  inherent 

hypothesis  that  is  always  made  when  dealing  with  images  modeled  as 

random  fields,  the  ergodicity  assumption.  This  in  turn  implies  that 

the  process  is  stationary  and  therefore  that  its  statistics  are 

invariant  with  translations  in  the  plane.  This  forces  some  non 

trivial  constraints  upon  the  different  nth  order  joint  density 

functions.  Without  going  into  two  much  details  (for  a  somewhat  more 

precise  description  see  [3]),  let  us  look  at  the  simple  case  of  the 

second-order  joint  density  functions  p. . (x. ,x.).  They  must  satisfy 

ID  i  D 


G-l 


S  Pii(Xi'Xn)  =  PC*-;) 

p  -*n  J  J 


x.=0 


(9) 


and 


G-l 

E 

v° 


pij  (xi'xj)  "  P(xi> 


(in) 


for  all  (i,j)  of  pixels  ( & j ) .  p(x)  is  the  first  order  density  of  the 
process  (the  same  at  every  point  because  of  the  stationarity) .  For 
higher  order  statistics,  more  complicated  relations  must  hold.  Since 
these  statistics  have  to  be  estimated  from  one  texture  sample  of 
finite  size  (using  the  ergodicity  assumption)  special  care  has  to  be 
taken  in  order  to  make  sure  that  constraints  such  as  Eqs.  (9)  and  (10) 
are  actually  verified. 
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Experimental  Results 

We  have  applied  our  algorithm  to  the  synthesis  of  two  natural 
texture  fields  "Water"  and  "Sand"  shown  in  Figure  2.  These  textures 
are  quantized  with  four  grey  levels  (G=4)  and  second  order  joint 
density  functions  (q=2)  are  estimated  in  a  15x15  neighborhood  (K=217). 
Results  of  the  synthesis  are  shown  in  Figure  3.  We  plan  to  experiment 
with  more  textures  and  different  values  of  the  parameter  G,  q  and  K  in 
the  very  near  future. 

Appendix  A 

We  want  to  minimize: 


&2  =  E{(xk~xk)2}  (A-l) 

with  respect  to  the  variables  where  x^  is  given  by  Eq .  (3). 

Setting  — 2.  equal  to  0  for  i=l,...,K  yields 


K 

E{E{3ck|yjL}(xk-2ajE{xk|yj})>  =  fi 


(A-2 ) 


and  therefore 


K 


2  a.E{x  |y.} 

j=l  3  *  3 


=  E{Xj^E{Xj,  |y^}  } 


(A-3 ) 


for  i=l,...,K. 

Appendix  B 

2 

Let  p  and  o  be  the  mean  and  variance  of  x^.  If  x^  is  given  by 
Eq.  (8),  3  is  adjusted  to  minimize 
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6'2  =  (xk-xk)2}-a2)2  (B-1) 

where  5  denotes  the  mean  of  It  is  readily  shown  that 

E{xk2}  =  B2+JTa5  =  g2+ST-S;  (R-2) 


and  therefore 


K 

EtT:  }  -  a. 

K  i=l  x 


(B-3) 


*2  = 


(g2+aT-b-M2  (£ 
i=l 


a,)2-a2)2 


(P-4  ) 


if  the  quantity  aT*£-y2  (.^ou)  2-o2  is  positive,  ^  is  minimized  for  g=D 
otherwise  for 


B2  =  a2+p2(Ea.)2-cT-S  (P-51 

i=l  1 
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2.8  Higher-Order  Texture  Synthesis  Models  and  Residual  Examination 
D.D.  Garber  and  A. A.  Sawchuk 


Introduction 


Earlier  work  [1]  has  shown  the  usefulness  of  the  Markov  and 
linear  autoregressive  model  in  the  analysis  and  synthesis  of  natural 
textures.  Each  model  was  assumed  to  be  the  random  process  used  to 
generate  a  texture.  Then,  based  on  this  assumption,  model  parameter 
estimates  were  made  using  data  from  a  parent  texture  and  these 
parameter  estimates  were  then  used  to  generate  simulations  of  that 
texture  using  the  assumed  model.  The  linear  model  was  developed  to 
add  simplicity  to  the  earlier  Markov  model  which  requires  a  large 
number  of  storage  locations  for  generation  parameters.  The 
simplification  arising  from  such  a  model  reduction  naturally  causes 
information  to  be  lost,  resulting  in  the  generation  of  less-appealing 
textures.  It  may  be  possible  to  regain  some  lost  information  by 
adding  important  parameters  to  the  reduced  model ,  thus  increasing  its 
complex ity. 

Markov  and  F i rst-Order  Models  ^ 

Texture  generation  work  up  to  the  present  Tl]  has  been  based  on 
the  concept  of  generating  a  pixel,  V  on  a  Hne  of  an  image  based 
on  the  values  of  previous  pixels  which  are  located  above*  or /to  the 
left  of  that  pixel  (Figure  1).  The  general  Markov-N  model  does  this 
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by  using  a  set  of  generation  parameters 


GV  *Vi'V2' — ,VN*  (1) 

N+l 


where 


GV  ' V2 ' *  *  * ^ ^N+l 

N+l 

and  P(A|B)  represents  the  conditional  probability  of  A  given  B.  Using 
this  model,  for  each  set  of  V.  ,  i=l,...,N  we  have  a  unique 
distribution  of  VN+1  which  can  assume  any  form.  Note  that  the  set  of 
parameters  (1)  not  only  determines  the  mean  and  standard  deviation  of 
the  next  pixel  VN+^  but  also  all  higher  moments.  The  limitation  of 
this  model  is  realized  rapidly  when  it  is  understood  that  gN+1  of 
these  parameters  must  be  used  to  generate  a  texture  as  each  unique  set 
of  v  ,  i=l,..,N+l,  requires  the  storage  of  a  generation  parameter 

where  g  is  the  number  of  distinct  grey  levels  a  pixel  can  assume  and  N 
is  the  number  of  pixels  in  the  generation  kernel  to  be  used  in  the 

also  known  as  an 

(2) 


this  case  each  VN+1 
depends  on  a  linear  combination  of  the  previous  and  some 
unexplained  variance  e.  In  this  case,  only  N  parameter  estimates  are 


generation  process. 


The  less-complex  linear  model  f 1 3 , 
autoregressive  model, 

Y  =  Xp+e 


where 


Y  =  V, 


N+l 


T 

X1  - 


N 
1 

was  then  proposed  to  generate  textures.  In 
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required.  (N  coordinate  points  are  also  required  if  the  kernel  is 
non-contiguous  [2]).  The  reduction  of  a  model  containing  g  estimates 
to  one  containing  only  N  will  naturally  cause  a  loss  of  information  to 
occur.  Geometrically  speaking,  we  are  approximating  an 
(N+l ) -dimensional  surface  with  an  N-dimensional  hyperplane.  In  the 
linear  model  case,  we  obtain  an  estimate  of  the  mean  of  VN+^  rather 
than  its  complete  distribution  for  a  given  set  of  V^,  i=l,...,N. 
Randomness  used  in  the  generating  process  involving  this  model  was 
obtained  by  assuming  that  £  has  some  fixed  distribution  with  zero 
mean.  The  extent  to  which  the  hyperplane  approximates  the 
(N+l ) -dimensional  surface  and  the  degree  to  which  the  eye  can  detect 
this  difference  will  determine  an  observer's  ability  to  discriminate 
between  a  texture  generated  by  the  Markov  model  versus  the  linear  one. 

The  degree  to  which  the  Markov  model  may  be  used  to  explain 
natural  texture  has  been  demonstrated  only  in  the  case  of  binary 
textures.  In  that  case,  the  model  worked  well,  exhibiting 
short-comings  mainly  in  cases  where  irregular  macro-structure, 
exceeding  the  size  of  the  neighborhood  of  V^s  used,  was  present. 
Application  of  this  model  to  continuous-tone  textures  remains  to  be 
undertaken  and  is  a  rather  large  task  in  itself  due  to  the  size  and 
amount  of  data  to  be  stored  and  utilized  in  the  form  of  generation 
parameters . 

Second-Order  Model 

When  we  say  that  a  model  is  linear  or  nonlinear,  we  are  referring 
to  linearity  or  nonlinearity  in  the  parameters.  The  value  of  the 
highest  power  of  an  independent  variable  in  the  model  is  called  the 
order  of  the  mod^l.  For  example, 

Y  =  61V1+6uv2+80+e  (3) 

is  a  second-order  linear  model.  A  general  second-order  linear  model 
with  two  independent  variables  may  be  written  as 


A  full  second-order  model  with  K  independent  variables  will  employ 
2 

(K  +3K)/2  terms  in  addition  to  the  3q (constant)  and  e  (error)  terms. 
Second-order  models  have  been  particularly  useful  in  studies  where 
surfaces  must  be  approximated  by  polynomials  of  low  order.  In  all 
cases,  a  second-order  model  will  "fit"  given  data  as  well  as  or  better 
than  a  first-order  model  containing  the  same  independent  variables  as 
the  first-order  model  is  a  subset  of  second-order  models.  This  does 
not  imply  that  the  second-order  model  will  be  more  correct  however  as 
the  process  which  we  are  attempting  to  explain  may  be  in  fact  a  linear 
first-order  process  or  some  other  type. 

It  is  also  important  to  note  that  the  covariance  of  the  are 

required  in  order  to  obtain  least-square  estimates  of  the  parameters 

3^  in  the  first-order  model  [2).  Covariance  is  essentially  a 

second-order  statistic  on  a  set  of  data.  Therefore,  estimating  the 

parameters  of  a  second-order  model  will  require  the  use  of 

fourth-order  statistics.  Specifically  the  correlation  of  terms 

and  will  be  utilized.  This  may  cause  serious  problems  as  in 

may  cases  the  variables  in  a  second-order  model  will  be  highly 
.  2 

intercor  related .  For  example,  the  terms  V^,  V-^  and  V-^V^  (if  is 

highly  related  to  V^)  may  be  strongly  correlated.  This  situation, 
often  referred  to  as  multicoil  inear ity,  may  cause  problems  during  the 
inversion  of  the  estimated  correlation  matrix,  a  necessary  step  in 
model  parameter  estimation.  For  this  reason,  care  should  be  exercised 
during  the  analysis  of  second-order  models. 

The  use  of  a  second-order  model  to  approximate  the  surface  of  the 
general  Markov  model  could  have  many  advantages  over  a  first-order 
model.  An  example  of  fitting  such  a  model  in  one  dimension  to  a  given 
set  of  data  is  shown  in  Figure  2. 


Still  the  linear  first-order  model  may  provide  a  good  fit  to  the 
data  and  the  magnitude  of  the  unexplained  variance  in  the  data  may  be 
large  enough  that  the  improvement  due  to  the  addition  of  second-order 
terms  to  the  model  may  be  barely  noticeable.  In  two  dimensions,  the 
fitting  problem  is  one  utilizing  a  quadric  surface  such  as  an  elliptic 
paraboloid  or  hyperbolic  paraboloid  versus  a  plane  to  fit  a  given  set 
of  data.  Again,  the  fit  may  or  may  not  be  markedly  better. 

An  Example 

In  order  to  test  whether  the  addition  of  second-order  terms  to 
the  linear  model  would  improve  the  quality  of  generated  textures,  the 
texture  sand  (Figure  3)  was  chosen.  It  was  felt  that  the  linear  model 
containing  only  first-order  terms  failed  to  accurately  simulate  this 
texture  (Figure  4)  and  that  improvement  could  be  made  using  a 
second-order  model. 

In  this  case  a  kernal  size  of  N*70  neighbors  was  chosen  based  on 

ideas  discussed  in  [2]  for  the  first-order  model  case.  That  is,  each 

pixel  generated  depended  on  the  linear  combination  of  70  pixels  above 

or  to  the  left  of  it.  When  adding  crossterms  to  the  model  it  is  most 

important  keep  in  mind  that  for  a  neighborhood  of  N  points,  /2 

second-order  terms  must  be  considered  for  possible  entry  into  the 

second-order  model.  For  N=70  this  is  2450  crossterms.  As  it  is 

numerically  infeasible  to  invert  a  ( 2450+70 ) x (2450+70 )  matrix  a  search 

procedure  known  as  stagewise  regression  was  used.  Basically,  one  uses 

the  best  first-order  linear  model  which  is  a  linear  function  of  M 

pixels  in  a  neighborhood  and  calls  that  the  new  variable  z^ .  Then  as 

many  crossterms  as  possible  are  entered  into  the  model  as  variables 

z .,  i=2,...,N  ,  where  N  is  a  more  manageable  figure  (less  than  100). 
i  s  s 

Model  parameter  estimates  are  obtained  at  this  step  and  the  best 
Nf<Ns,  (those  chosen  for  entry  into  the  model  according  to  some 
criterion  [2])  are  kept  in  the  model  while  the  remaining  NS_N|  are 
discarded.  At  the  next  step,  N  -N  ,  crossterms  which  have  not  been 

ST 

previously  tested  are  examined  for  entry  into  the  model.  Again,  the 
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best  N  are  chosen  and  the  rest  are  discarded.  The  process  continues 
until  all  possible  crossterms  of  interest  have  been  tested  for  entry 
into  the  model.  Keeping  N  f  large  helps  to  insure  that  valuable  model 
terms  are  not  discarded  during  the  process.  Keeping  N  small, 
however,  reduces  the  number  of  steps  to  be  performed,  as  it  allows  the 
examination  of  more  new  terms  at  each  step.  Each  step  involves  the 
computation  and  inversion  of  a  large  covariance  or  correlation  matrix. 
Thus  there  is  a  trade-off  to  be  made.  This  procedure  does  not  provide 
a  true  least-squares  solution  for  the  variables  and  cross-terms 
entered  into  the  final  model.  Still  the  solution  is  a  very  good  one 
in  most  cases  and  is  probably  the  best  available  given  the 
dimensionality  constraints  of  our  problem. 

An  analysis  of  this  type  was  performed  on  the  sand  texture  for 
N=7 0 .  All  cross-products  containing  these  points  were  investigated. 
The  results  of  the  synthesis  using  this  second-order  model 
incorporating  white  gaussian  noise  with  mean  zero  and  fixed  variance 
to  simulate  the  error  term  of  the  model,  e,  are  shown  in  Figure  5. 
The  results  show  only  a  slight,  almost  undetectable  improvement-much 
less  than  desired.  This  would  indicate  that  the  improvement  of 
texture  synthesis  due  to  the  addition  of  second-order  terms  into  the 
linear  model  is  minimal  in  this  case. 

Examination  of  the  Residual  Image 

A  residual  image  is  that  image  which  is  formed  by  the  differences 

A  A 

Vn+1-Vn+1  where  vn+1  is  an  observed  pixel  value  and  VN+1  is  the 
corresponding  fitted  value  obtained  by  use  of  the  linear  model  (?) . 
We  can  see  that  these  residuals  are  the  differences  between  what  is 
actually  observed  and  what  is  predicted  by  the  model,  that  is,  the 
amount  of  variation  which  the  linear  model  has  not  been  able  to 
explain.  The  usual  assumption  about  the  errors,  e,  in  our  model  (2) 
is  that  they  are  independent  with  zero  mean  and  constant  variance.  An 
additional  assumption  of  normality  allows  us  to  make  F-tests  when 
checking  variables  and  crossterms  for  significance  [21.  If  our  fitted 
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model  is  correct,  the  residuals  should  exhibit  tendencies  that  confirm 
these  assumptions  or  at  least  should  not  exhibit  a  denial  of  the 
assumptions . 

A  time  or  space  sequence  plot  of  residuals  5s  most  fitting  in  the 
case  of  two-dimensional  texture  synthesis.  This  is  essentia] ly  an 
image  containing  the  pixel  differences  or  residuals  vn+i-vn+1* 
Naturally  we  would  expect  these  errors  to  be  small  as  merely 
subtracting  one  pixel  from  its  nearest  neighbor  would  yield  a  smal] 
value  in  most  natural,  low-noise  images.  Such  an  image  of  residuals 
was  generated  for  the  sand  and  linearly  rescaled  to  show  the  detai] 
present  in  the  image  (Figure  6).  Definite  patterns  are  seen  to  exist 
in  this  image  and  thus  a  violation  of  the  independence  assumption  is 
indicated.  Ideally,  this  residual  image  would  be  uncorrelated  noise. 

A 

A  histogram  showing  the  number  of  V  ,  occurring  at  each  pixel 

^  N+l 

value  is  shown  in  figure  7.  A  plot  indicating  the  mean  of  the 

A 

residuals  V  -V  „  versus  V  is  shown  in  figure  8.  As  would  be 
N+l  N+l  N+l 

expected,  residuals  where  the  V  .  is  less  than  0  will  have  a  mean 

N+l 

less  than  zero  and  those  residuals  where  the  V  ,  is  greater  than  255 

N+l 

will  be  likewise  positive.  Figure  9  shows  a  similar  plot  of  the 

A 

standard  deviation  of  the  residuals  versus  V  .  These  three  figures 

N+l 

seem  to  indicate  that  the  distribution  of  the  error  in  the  model  is 

A 

related  to  the  value  V  Therefore  the  assumption  of  constant  error 

variance  is  questionable.  It  may  be  reasonable  to  drive  the 
generation  process  with  noise  which  does  not  have  stationary  mean  or 
variance.  The  effect  of  such  a  change  in  the  generation  process 
remains  to  be  studied.  It  is  also  not  yet  clear  what  effect  such  a 
change  will  have  on  a  texture  synthesis.  Figure  10  shows  the 

A 

distribution  of  error  V  -V  ,  (a  histogram  of  the  residual  image). 

Given  the  apparent  violation  of  assumptions  concerning  the  model 
we  must  re-examine  whether  these  violations  may  be  explained  in  some 
specific  way  and  whether  they  can  be  adjusted  for  in  a  proper  manner. 
Such  adjustments  may  involve  transformations  of  existing  variables, 
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the  generation  of  new  theory  concerning  texture  synthesis  or  the  use 
of  non-stat iona ry  non-random  noise  in  the  generation  process.  These 
ideas  and  others  remain  to  be  investigated. 

Conclusions 


More  investigation  into  the  effect  of  the  addition  of 
second-order  terms  into  the  texture  model  must  be  done.  More 
conclusive  results  will  only  be  developed  after  this  study  has  been 
performed  on  a  variety  of  textures.  The  study  of  the  residual  image 
indicates  that  many  important  model  assumptions  are  being  violated  and 
thus,  at  least  in  some  ways,  our  model  may  be  inadequate. 
Particularly,  in  the  case  of  sand,  the  texture  generation  method  fails 
to  generate  irregularly  shaped  sharp  edges,  therefore  it  may  be 
beneficial  to  propose  changes  into  our  model  which  will  induce  these 
edges . 
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2.9  Computer  Analysis  of  Moving  Images 
B.  Bhanu  and  O.D.  Faugeras 


Introduction 


In  the  past  image  understanding  systems  capable  of  doing  simple 
image  description  tasks  for  a  narrow  class  of  images  have  been 
constructed.  The  development  of  systems  for  the  purposes  of 
detections,  identification,  localization  and  tracking  the  objects 
through  a  sequence  of  images  will  be  a  very  desirable  and  important 
features  of  such  machine  vision  systems.  For  example,  in  the  sequence 
the  objects  of  interest  may  be  the  enemy  targets  tracked  by  the 
imaging  missile  seekers  which  incorporate  sophisticated  signal 
processing  to  track  the  targets,  aircrafts  or  blood  cells  etc.  The 
constituent  of  the  scene  may  be  moving  in  two  or  three  dimensional 
space  against  a  stationary  foreground/background  and  sequence  of 
images  may  be  in  the  optical  or  some  other  band  of  the  spectrum.  The 
task  at  hand  is  to  consider  the  incorporation  of  various  diverse 
sources  of  knowledge  that  lead  to  the  science  of  image  understanding. 
Recently  both  temporal  and  spatial  sequences  of  images  have  been  used 
in  understanding  and  analyzing  a  wide  variety  of  processes  in  physics, 
biology,  medicine,  meterology,  navigation  etc.  11-5].  Special 
purpose  systems  have  been  built  to  record  and  store  the  sequences  of 
television  frames  and  to  correlate  information  from  two  cameras  as 
well  as  from  a  sequence  of  images  [5,6]. 

Most  of  the  work  in  the  field  of  image  processing  has  been 
concerned  with  a  single  frame,  where  the  primary  objective  is  to 
develop  methods  to  detect,  represent,  store  and  manipulate  the  spatial 
features  of  a  scene.  Although  a  sequence  of  images  is  basically  a 
sequence  of  static  images  with  a  given  time  function  relating  the 
order  and  time  interval  between  the  images  of  the  sequence,  yet  in 
analyzing  the  sequence  of  images,  not  only  information  must  be 
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extracted  from  each  frame,  but  also  from  the  sequence  as  such  so  that 
a  description  of  the  sequence  can  be  obtained.  Noise,  occlusion, 
perspective  changes  of  objects,  appearance  of  new  objects, 
disappearance  of  old  objects,  shape  changes  further  complicate  the 
tracking  and  matching  problem. 

In  this  paper  we  consider  a  new  dynamic  image  understanding 
system  for  analyzing  the  behavior  of  2-D  curvilinear  objects  through  a 
temporal  sequence  of  images.  Earlier  model  based  systems  are 
principally  hierarchical  image  understanding  systems  dealing  with 
polygonal  shaped  objects  [8]  or  curvilinear  objects  [2,9]  or  the 
objects  from  blocks  world.  They  suffer  a  lack  of  adaptability  to  the 
task  at  hand,  require  more  processing  and  consider  only  few  objects 
whose  number  remain  fixed  in  the  sequence.  The  larger  number  of 
objects  exhibit  more  complex  interactions  among  the  objects  and  the 
disappearance  and  appearance  of  objects  through  the  boundaries  of  the 
images  enhances  the  difficulty  of  the  problem  and  require  the  control 
on  symbols  and  features  of  the  objects  in  a  feedback  loop. 

The  next  section  describes  in  brief  the  previous  work,  followed 
by  a  description  of  the  system  for  recording  the  sequences  and  types 
of  scenes.  Finally,  we  consider  the  dynamic  image  understanding  model 
for  analyzing  the  motion  and  shape  of  the  constituent  of  the  images. 

P revious  Work 

At  present  digital  as  well  as  optical  image  processing  methods 
are  being  used  for  the  motion  analysis  purposes.  Both  signal  and 
symbolic  approaches  are  under  investigation.  Originally,  the  work  was 
stimulated  by  the  interest  of  analyzing  the  movement  of  clouds  from  a 
sequence  of  satellite  images.  Work  has  proceeded  from  there  to  the 
analysis  of  abstract  patterns  in  two  dimensions  and  the  analysis  of 
scenes  from  the  blocks  world  [2,3,5].  Variety  of  techniques  have  been 
used  by  the  researchers  working  on  the  interframe  coding  for  comparing 
image  frames,  segmenting  them  into  stationary  and  non-stat iona ry  areas 
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and  estimating  the  translational  velocity  of  moving  object.  A  number 
of  similar  techniques  have  also  been  used  by  the  researchers  in  the 
medicine,  biology,  industrial  automation,  robotics,  and  other  related 
areas.  Research  efforts  involved  in  the  analysis  techniques  for  image 
sequences  can  be  classified  into  the  following  categories. 

A.  Techniques  based  on  cross-correlation 

B.  Centroid  matching 

C.  Change  detection  and  Image  differencing 

D.  Dynamics  scene  modeling 

E.  Velocity  Estimation  techniques 

F.  Hough  Transform  techniques 

G.  Optical  methods  and  tracking 

H.  Moving  target  indication  filter 

I.  Kalman  filtering  techniques 

These  techniques  have  been  discussed  in  detail  in  [2,3,7]. 

System  Configuration  and  Input  Scenes 
1 .  System  Configuration: 


For  acquiring  both  real  and  synthetic  sequences  of  images,  a  real 
time  video  acquisition  and  digital  display  System  shown  in  Fig.  1  has 
been  used.  The  system  was  built  at  the  USC  Image  Processing  Institute 
recently  [6].  The  system  under  the  control  of  the  PDP  11/34  (RT-11 
Operating  system)  processor  will  acquire,  digitize  and  store  a  number 
of  consecutive  subframes  of  video  images  from  the  monochrome  TV 
camera.  During  and  after  the  acquisition  phase,  the  acquired  image 
can  be  displayed  on  the  monochrome  TV  monitor.  After  the  acquisition 
phase,  the  acquired  image  subframes  can  be  transferred  from  the  image 
memory  (512x512x8  bits)  to  the  PDP  11/34.  An  image  or  image  subframes 
can  also  be  transferred  from  the  PDP  11/34  to  the  image  memory.  The 
maximum  size  of  the  image  subframe  (window  area)  is  140x400  pixels. 
The  minimum  sampling  interval  between  two  subframes  is  of  the  order  of 


Monitor 


Figure  1.  Real  time  video  acquisition  and  digital  display  system 


0.75  seconds.  It  is  limited  by  the  size  of  the  subframe  and  the 
amount  of  time  required  to  transfer  and  store  the  image  into  a  file  on 
the  disk  or  tape. 

Programs  have  been  written  by  authors  on  PDP  11/34  for  the  system 
shown  in  Fig.  1,  which  allow  the  user  to  store  the  subframes  of 
specific  sizes  at  a  specified  rate.  Odd  and  even  lines  in  the 
subframe  can  be  stored  under  separate  files  on  the  disk  for  the 
interlacing  to  be  carried  out  later  on  or  they  can  be  interlaced 
before  storing  it  on  the  disk.  Images  so  stored  can  be  put  back  on 
the  display  under  the  control  of  PDP  11/34  for  observing  a  part  of  the 
sequence.  Sequences  of  the  images  obtained  are  transferred  on  PDP-1 0 
(Tenex  operating  system)  via  a  format  conversion  from  RT-11  to  RSX-11 
and  further  processing  is  done  on  Tenex. 

2 .  Input  Scenes 

Using  the  system  described  above  we  have  generated  a  number  of 
synthetic  and  real  sequences. 


Synthetic  sequences  exhibit  the  motion  of  planar  curvilinear 
two-dimensional  objects  against  the  different  types  of  background. 
The  complexity  of  the  scene  in  these  sequences  vary  from  8  to  15 
objects.  For  generating  these  sequences  the  following  considerations 
have  been  taken  into  account. 

1)  Linear  and  rotational  motion  of  the  objects  is  allowed.  Some 
of  the  objects  may  not  move  at  all,  while  some  others  may 
show  random  behavior. 

2)  No  more  than  two  to  four  objects  may  occlude  each  other  at 
anytime  and  no  more  than  two  to  three  objects  may  appear  or 
disappear  through  the  boundaries.  At  the  most  an  object  may 
divide  itself  into  two  objects.  Thus  the  number  of  objects 
may  not  remain  fixed  throughout  the  sequence. 

3)  Objects  don't  have  holes  and  slow  changes  in  the  shape  of 
objects  and  background  are  allowed. 


4)  Changes  are  slow  from  frame  to  frame.  They  are  not  drastic 
or  discontinuous.  Of  course,  significant  changes  may  or  may 
not  be  present  between  the  two  frames  in  the  sequence. 

Figure  2  shows  some  frames  from  two  such  sequences  with  different 
backgrounds . 

In  order  to  carry  out  the  motion  and  shape  analysis  on  some  real 
sequences  of  images,  we  have  also  carried  out  experiments  at  the  USC 
medical  school,  to  analyze  the  behavior  of  lymphocytes  in  the  presence 
of  cancer  cells  in  vitro.  The  TV  camera  in  Fig.  1  is  connected  to  a 
NIKON  phase  contrast  microscope.  A  confluent  monolayer  of  human  skin 
cancer  cells  was  grown  on  microslides.  Using  this  slide  a  Sykes-Moore 
chamber  was  formed  and  filled  with  a  RPMI  1640  media  with  10%  serum. 
Lymphocytes  were  isolated  from  fresh  human  blood  by  gelatine  technique 
and  their  concentration  was  about  90%.  A  small  number  of  these 
lymphocytes  were  injected  into  the  chamber.  The  whole  system  was 
incubated  at  37°c.  The  movement  of  the  lymphocytes  was  recorded  at  30 
sec.  to  60  sec.  intervals  in  the  images  (128x128  pixels)  over  a 
period  of  about  24  hours.  128x128  pixel  image  corresponds  to  an  area 
of  60  ym  x  60  ym  of  the  slide  at  x200  magnification.  Some 

representative  pictures  are  shown  in  Fig.  3. 

Dynamic  I  mage  Understand i ng  Model 

For  the  motion  and  shape  analysis  of  objects  in  an  adaptive 
manner  we  consider  a  new  dynamic  image  understanding  model  presently 
under  investigation.  Its  simplified  block  diagram  is  shown  in  Fig.  4. 

As  seen  from  Fig.  4  that  this  model  achieves  the  symbols, 
features  and  segmentation  control  by  distributed  feedback.  It  is 
similar  to  the  heterarchical  model  suggested  by  Minsky  &  Papert  and 
used  by  Shirai  [10]  for  obtaining  line  drawings.  On  the  basis  of  the 
fact  that  the  predicted  model  of  the  scene  at  a  particular  time 
matches  with  the  scene  at  that  time,  we  do  simple  or  complicated 
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Fig.  3.  A  real  sequence  of  16  images. 
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processing  for  motion  as  well  as  shape  analysis, 
varies  with  the  dyramics  of  the  scene. 


Also  segmentation 


Throughout  this  study  emphasis  will  be  given  on  the  iterative, 
parallel  processes,  known  as  "relaxation"  techniques  [11-13].  In 
general,  "relaxation"  is  an  iterative  approach  for  classifying  a  set 
of  interrelated  objects.  In  "discrete  relaxation"  a  set  of  possible 
class  names  is  initially  associated  with  each  object.  At  subsequent 
Iterations,  class  names  are  discarded  from  an  object  if  they  are 
inconsistent  with  the  surviving  class  name  possibilities  for  other, 
related  objects;  this  is  done  for  all  objects  simultaneously.  The 
process  is  repeated  until  no  further  changes  can  take  place.  In 
"probabilistic  relaxation",  a  set  of  estimates  of  class  assignment 
probabilities  is  initially  associated  with  each  object.  At  subsequent 
iterations,  the  probabilities  are  adjusted  in  accordance  with  the 
support  they  receive  from  the  class  probabilities  of  related  objects. 
This  process,  when  repeated  often  leads  to  a  marked  reduction  in  the 
ambiguity  of  classification. 

Ullman  [10]  and  Faugeras  and  Berthod  [62]  consider  the  relaxation 
processes  as  problems  of  computing  constrained  optimization  by  local 
processes.  The  algorithm  of  Rosenfeld  et.  a  1 .  T 1 2 ]  uses  no  explicit 
measure  of  ambiguity  of  classification.  Faugeras  and  Berthod 
formulate  the  relaxation  algorithm  in  terms  of  minimizing  a  criterion 
function  which  includes  the  inconsistency  and  ambiguity  of 
classi f ication . 

Barnard  and  Thompson  [14]  and  Ullman  [10]  consider  the 
appl  ication  of  relaxation  algorithms  in  the  matching  of  two  images. 
In  [14]  an  initial  network  of  possible  matches  between  the  two  sets  of 
candidates  in  two  images  is  constructed.  Each  possible  match 
specifies  a  possible  disparity  of  a  candidate  point  in  a  related 
reference  image.  An  initial  estimate  of  the  probability  of  each 
possible  disparity  is  made,  based  on  the  similarity  of  subimages 
surrounding  the  points.  These  estimates  are  iteratively  improved  by 
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using  the  relaxation  algorithm  similar  to  of  Rosenfeld  et.  al.  [12]. 


There  are  a  number  of  drawbacks  in  merely  using  the  matching 
approach  for  analyzing  the  motion  of  objects  in  a  sequence  of  images. 
First,  although  such  matching  approaches  incorporate  spatial 
correlation  within  a  image,  they  do  not  consider  temporal  correlation 
at  all.  Since  the  successive  images  are  highly  correlated,  the 
inclusion  of  this  information  should  be  helpful  in  "better"  and  "fast" 
matching.  Second,  if  such  a  matching  approach  is  applied-  to  a  large 
number  of  frames,  it  will  be  very  expensive  computationally. 

Study  removing  such  drawbacks  and  the  details  of  the  model  along 
with  the  computational  issues  related  to  analyzing  the  sequences  of 
images  are  being  worked  out  and  will  be  reported  in  the  future. 
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2.10  Segmentation  of  Images  having  Unimodal  Gray  Level  Distributions 
B.  Bhanu  and  O.D.  Faugeras 


Introduction 


Image  segmentation  is  a  transformation  from  a  point  by  point 
representation  of  an  image  to  a  representation  of  the  image  as  a 
collection  of  meaningful  regions.  In  image  analysis  it  is  the  most 
primitive  process,  because  a  segmented  image  requires  less  storage 
than  the  original  and  can  be  analyzed  on  the  basis  of  shapes,  features 
and  other  characteristics  of  its  regions.  Various  approaches  based  on 
thresholding  have  been  used  by  many  researchers  for  the  segmentation 
of  both  monochrome  and  color  pictures  Tl-4].  A  good  survey  of  global, 
local  and  dynamic  threshold  selection  methods  is  given  in  fll. 

Normally,  in  the  application  of  these  techniques,  the  gray  level 
histogram  of  the  image  shows  two  or  more  peaks  and  the  preprocessing 
is  done  to  improve  the  histogram,  if  needed.  Local  properties  are 
also  used  to  compute  the  global,  local  or  dynamic  thresholds. 
However,  if  the  gray  level  histogram  of  the  image  is  unimodal,  then 
the  application  of  such  methods  give  poor  segmentation  results  and 
there  is  no  criteria  for  automatic  threshold  selections.  In  T 8 , 9 1 
Rosenfeld  and  Davis  discuss  the  histogram  modification  iteratively. 
Various  pictures  shown  in  [8]  already  consist  of  2  or  more  peaks.  Tn 
this  paper  we  consider  the  application  of  relaxation  methods  for 
segmenting  monochrome  pictures  having  unimodal  gray  level 
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distributions.  It  is  found  that  the  gradient  relaxation  method  gives 
very  good  segmentation  and  provides  the  automatic  selection  of 
thresholds.  Results  are  illustrated  with  the  aid  of  examples. 

Segmentation  Schemes 

Figure  1  shows  two  128x128  pixels  8-bit  images.  Image  in 
fig.  1(a)  has  been  obtained  using  the  experimental  setup  discussed  in 
another  paper  on  "Computer  Analysis  of  Moving  Images"  of  this  report. 
The  background  in  this  image  consists  of  a  confluent  monolayer  of 
human  skin  cancer  cells  and  the  small  circular  shaped  objects  are 
human  lymphocytes  and  red  blood  cells.  The  objective  is  to  get  the 
boundaries  of  all  the  cells.  Image  in  fig.  1(b)  is  taken  from 
DMA3.2048.  Here  the  objective  is  to  detect  roads  etc.  Gray  level 
histogram  of  these  images  are  shown  in  Fig.  2.  Note  that  the 
histograms  are  Gaussian  shaped.  Since  the  histograms  are  unimodal , 
there  is  no  automatic  method  for  segmenting  these  pictures. 

Commonly  used  difference  operations  such  as  gradient,  Laplacian 
and  Sobel  were  applied  on  the  images  shown  in  Fig.  1.  We  also 
considered  the  methods  based  on  thresholding  the  histogram  of  the 
picture  where  the  gradient,  Laplacian  and  edge  values  are  high  Til. 
However,  the  picture  so  obtained  has  unimodal  histogram  and  lacks  the 
criterion  for  segmenting  it  at  the  valley  of  two  peaks.  Thresholding 
at  the  gray  level  corresponding  to  the  mode  or  mean  of  the  filtered 
histogram  gave  very  poor  results.  Edge  detection  has  also  been  done 
by  convolving  the  image  in  Fig.  1  with  5x5  masks  corresponding  to  the 
ideal  step  edges  in  six  directions  as  described  in  [5],  The  magnitude 
of  the  convolved  output  and  the  direction  of  the  mask  giving  the 
highest  output  at  each  pixel  are  recorded  as  edge  data.  The  magnitude 
image  does  not  show  good  segmentation.  A  number  of  bar  masks  of 
various  sizes  and  orientations  have  also  been  used,  but  the  results 
were  poor  (for  the  image  shown  in  Fig.  la  the  width  of  an  edge  is 
about  5-6  pixels) . 
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Two  typical  128x128,  8  bit  images. 

(a)  Cell  image,  (b)  Image  taken  from 
DMA  3.2048 


Gray  level  histogram  of  the  images  in  Fig 

(a)  Histogram  of  the  image  in  fig.  1(a). 

(b)  Histogram  of  the  image  in  fig.  1(b). 


Mean=163. 76 
Mean=26 . 47 . 
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Segmentation  Using  Relaxat ion  Methods 

Having  considered  the  basic  differencing,  edge,  bar  operators  and 
their  variations,  we  focussed  our  attention  on  the  probabilistic 
relaxation  methods  for  segmentation  [6,7],  "Relaxation"  is  an 
iterative  approach  to  classifying  a  set  of  interrelated  objects.  In 
"probabilistic  relaxation,"  a  set  of  estimates  of  class  assignment 
probabilities  is  initially  associated  with  each  object.  At  subsequent 
iterations,  the  probabilities  are  adjusted  in  accordance  with  the 
support  they  receive  from  the  class  probabilities  of  related  objects. 
This  process,  when  repeated  often  leads  to  a  marked  reduction  in  the 
ambiguity  of  classification. 

Nonlinear  Probabilistic  Relaxation  Algorithm 

First  we  consider  the  nonlinear  probabilistic  relaxation 
algorithm  proposed  in  [6]  which  has  been  extensively  used  in  edge 
detection,  edge  and  line  enhancement,  clustering  etc.  In  order  to  use 
this  algorithm  for  segmentation  purposes,  suppose  we  have  a  set  of  N 
objects  (pixels)  ap,..^,  which  fall  into  two  classes  X and  X2 
corresponding  to  the  white,  and  black  gray  values.  For  the  same 
neighborhood  of  the  two  classes  (known  as  homogeneous  case)  ,  the 
probability  of  being  white  at  the  (n+l)th  iteration  is  given  by. 


Pf+1) 


<V 


Pjn)(V^ni  <V 


E  p in)  i v-sf’ < V 


k  =  1,2 

i  =  1, . . .  ,N 


(1) 


where  the  consistency  vector 


q. 


is, 
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a^j's  are  called  as  compatibility  coefficients.  Considering  the 
neighborhood  \T  of  8  around  the  central  pixel  a^for  both  classes  and 
taking  equal  values  of  ct^'s  Eq.  (2)  becomes. 
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Pij(Xk|X^)  is  the  conditional  probability  that  object  a^  belongs  to 
class  Xk  given  that  object  a^  e  vi  belongs  to  class  X£.  For  the  two 
class  problem,  we  take 

Pi-jUilV  -  Piju2|x2)  -  i  ,4, 


Thus , 


pij^A1^2)  =  pij(*2^1*  =  0 


If  ^ Al^  8  pi^Al^ 
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Noting  that 
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Eq.  (1)  can  be  written  as. 
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So  at  each  iteration  given  the  values  of  p^fA^),  we  compute  qfm)(A^) 
from  Eq.  (6)  and  then  Pi(m+1)(x1)  from  Eq.  (7). 

In  the  case  when  we  consider  the  neighborhood  of  8  for  class  X1 
and  the  neighborhood  4  for  class  X]/  Eq.  (6)  becomes. 
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(8) 


The  initial  assignment  of  probabilities  (X,)  can  be  done  by  simply 
normalizing  the  gray  values  of  the  image  by  255.  However,  it  does  not 
take  into  account  the  mean  value  of  the  picture.  So  the  following 
transformation  is  used, 


FACT  * 


,1 -Global 
255 


mean 


)  +  0.5 


(9) 


where  I  =  intensity  value  at  a  pixel  and  FACT  is  a  constant  whose 
value  is  determined  as  follows: 


Fact  =  1 


if  I  >  Global  mean 


(10) 


=0.7  to  0.9  if  I  <  Global  mean 


In  Eq.  (9)  if  Fact*  (l^i°b|l_mean }  <_0^5f  then  p^O)^)  is  taken  as 
zero  at  that  pixel. 


This  transformation  is  motivated  by  the  fact  that  the  gray  level 
histogram  of  a  typical  natural  image  is  usually  skewed  toward  the 
darker  gray  values  and  therefore,  the  majority  of  the  pixels  possess 
gray  value  less  than  average.  Figs.  3  and  4  illustrate  the 
application  of  this  method  to  the  images  shown  in  Fig.  1.  As  can  be 
seen  from  the  histogram  plots,  in  the  first  few  iterations  a  sort  of 
histogram  equalization  takes  place  and  then  black  and  white  clashes 
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(a)  Iteration  1 


(c)  Iteration  3 


(e)  Iteration  4 


Flg*  3.  Results  of  Nonl 
iterations  and 
in  fig.  l (a) . 


(i)  Iteration  7 


(j)  Histogram  of  fig.  3 ( i ) 


(k)  Iteration  9 


(1)  Histogram  of  fig.  3 (k) 


Fig.  3.  (CONT.) 


Histogram  of  fig.  4 (a) 


Results  of  Nonlinear  relaxation  method  at  various 
iterations  and  corresponding  histograms  for  the  image 
in  fig.  1(b).  FACT=1  (see  text). 
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get  separated.  However,  in  this  method  convergence  is  slow, 
uncontrolled  and  does  not  lead  to  the  automatic  selections  of 
threshold  at  the  valley  of  two  peaks  which  are  not  far  apart.  Results 
obtained  by  thresholding  at  the  mean  are  discussed  in  the  next 
section.  Gradient  relaxation  method  discussed  next  overcomes  these 
disadvantags. 

Gradient  Relaxation  Algorithm 

Faugeras  and  Berthod  [7]  proposed  a  different  approaches  for 
relaxation  algorithm  which  is  based  upon  the  explicit  use  of 
consistency  and  ambiguity  to  define  a  global  criterion  upon  the  set  of 
objects.  The  criterion  is  minimized  using  the  projection  gradient 
techniques.  In  [7]  they  consider  the  minimization  of  the  criterion, 

C3  =  aC1+(l-a)C2  (11) 

where,  a  varies  between  0  and  1,  and  weights  the  relative  importance 
we  assess  to  consistency  versus  ambiguity,  and 

N 

C1  =  2N  H  2  (12) 
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the  present  study  we 
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and  carry  out  its  maximization  using  the  projection  gradient  approach. 
This  criterion  is  easier  to  manipulate  computationally. 

For  the  two  class  case  A^,  A2  and  the  same  8-neighborhood  for 
these  classes  and  assuming  the  same  conditional  probabilities  as  in 
the  earlier  discussion  on  nonlinear  relaxation  method,  gradient  of  the 
criterion  C  is  obtained  as, 


3Pi(Ai)  "  2cIi(xi> 


■5---  -,V  v  =  2q.  (A_) 
3pi(A2)  2 


The  initial  assignment  of  probabilities  is  taken  to  be  the  same 
as  in  the  nonlinear  relaxation  method.  The  iteration  of  p^'s  is  given 


(n+1)  . ,  .  _  (n)  .  p(n)r  dC' 1 
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(18) 
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Pi  (X2)  -  Pi  (x2)+PiPi 


In  order  to  have  p x  ^ +p ^(  x 2)  =  1»  the  projection  of  the  gradient 
should  be  such  that  [7], 
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and 
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So  the  iteration  eqs .  in  (17)  and  (18)  reduce  to 

pl"+llUl>  -  P^’up+P^uqilxp-I)  (21! 

p^n+l) (X^j  _  (i2)+p^n) [2qi (X2)-1J  (22) 

Normally,  p|n^  is  kept  constant  for  all  pixels  during  each  iteration 
and  is  determined  to  have  the  largest  possible  value  such  that  p^ * s  at 

(n+l)th  iteration  still  lie  in  the  bounded  convex  region  of  2N 

dimensional  Euclidean  space.  However,  in  the  2  class  considered  it  is 

easier  to  compute  for  each  pixel  (although  we  may  not  be 


folio*  ing  the  gradient  exactly)  and  is  obtained  from  Eq .  (21)  and  (2?) 
as , 
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where  K=l,2,  and  and  are  constants  less  than  1.  These  values 

may  be  selected  to  bias  a  class  and  control  the  speed  of  convergence. 
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Figs.  5  and  6  show  results  obtained  at  various  iterations  and 
corresponding  histograms.  Note  that  at  the  first  iteration  itself  we 
get  two  peaks  separated  by  a  valley  which  can  be  used  to  automat i ca 1 1 y 
select  the  threshold  to  obtain  segmentation.  As  the  number  of 
iteration  increases  the  two  peaks  get  apart,  average  brightness 
increases,  and  the  convergence  of  probabilities  takes  place  as 
expected.  When  the  peaks  are  far  apart,  thresholding  can  be  done  at 
the  mean  value.  Segmentation  result  obtained  by  semi-thresholding  the 
gray  level  histogram  are  shown  in  Figs.  7  and  8.  These  clearly  show 
that  gradient  relaxation  method  gives  better  segmentation  than 
nonlinear  relaxation.  For  the  nonlinear  relaxation  case,  we  have 
shown  the  segmentation  results  only  when  the  two  peaks  are  separated. 

In  the  above  example  we  have  considered  the  8-ne ighborhood  for 
both  classes  A ^  and  x^.  If  we  consider  the  8 -ne ighborhood  for  class 

\  and  4-neighborhood  for  class  A-,,  the  gradients  are  obtained  as. 
i  2- 


(b)  Histogram  of  fig.  5 (a 


(e)  Iteration  4 


(f)  Histogram  of  fig.  5(e) 


Fig.  5.  Results  of  Gradient  Relaxation  method  at  various  iterations 
and  corresponding  histograms  for  the  image  in  fig.  1(a). 
FACT=0 . 9  ,  (*2=0.2,  a 2=0.1  (see  text). 


(b)  Histogram  of  fig.  6(a) 


(d)  Histogram  of  fig.  6(c) 


(f)  Histogram  of  fig.  6(e) 


Fig.  6.  Results  of  Gradient  Relaxation  method  at  various  iterations 
and  corresponding  histograms  for  the  image  in  fig.  1(b). 
FACT=1,  a^=0.1,  (*2  =  0.05  (see  text). 
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(a)  Semi-thresholding  of  fig.  5(a)  at 
the  valley  bottom  (see  fig.  5(b)) 


(b)  Semi-thresholding  of  fig.  5(g) 
at  the  mean 


(c)  Semi-thresholding  of  fig.  3(g) 
at  the  mean 


(d)  Semi-thresholding  of  fig.  5 ( i ) 
at  the  mean 


(e)  Semi-thresholding  of  fig.  3 ( i ) 
at  the  mean 


Fig.  7.  Segmentation  results  for  Gradient  and  Nonlinear  relaxation  methods 
obtained  by  semi-thresholding  at  various  iteration  for  the  image 
in  fig.  1(a). 


(a)  Semi-thresholding  of  fig.  6(a) 
at  the  valley  bottom  (see 

fig-  6(b)) 


(c)  Semi-thresholding  of  fig.  6(e) 
at  the  mean 


(d)  Semi-thresholdinn  of  fig. 
at  the  mean 


Fig.  8.  Segmentation  results  for  Gradient  and  Nonlinear 

relaxation  methods  obtained  by  semi-thresholdinq  at 
various  iterations  for  the  image  in  fig.  1(b). 
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Now  the  projection  of  the  gradient  can  be  obtained  as  before. 
The  application  of  the  iteration  thus  obtained  gave  the  results  very 
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much  like  those  obtained  using  Eq.  (21). 

This  method  has  been  applied  on  a  number  of  aerial  and  other 
pictures  having  unimodal  distributions  and  gave  good  segmentation 
results . 

Conclusions 


It  has  been  shown  that  the  gradient  method  of  relaxation  can  be 
used  as  a  tool  for  segmentation  of  black  and  white  pictures  having 
unimodal  distributions.  This  method  provides  automatic  selection  of 
the  threshold  which  leads  to  good  segmentation  results. 

In  this  study  we  have  considered  only  two  classes.  Tn  scenes 
with  many  different  objects  as  in  the  case  with  aerial  photographs, 
there  are  more  classes  although  the  histogram  may  have  only  one  peak 
because  the  range  of  intensities  for  each  object  will  probably  overlap 
with  the  ranges  of  other  objects.  Two  generalizations  are  possible. 
In  the  first  we  simply  extend  the  method  for  more  classes,  which  are 
to  be  known  a  priori.  In  the  second  case,  we  apply  the  procedure 
outlined  in  this  paper  to  each  segmented  region  iteratively  in  a  tree 
structure  until  the  area  of  a  region  has  reduced  to  a  specified  limit. 
If  a  region  is  homogeneous  in  intensity,  then  we  will  not  expect  to 
find  two  peaks  and  will  stop  iterating  after  a  predefined  number  of 
iterations.  An  experiment  investigating  the  sensitivity  of  this 
method  to  homogeneous  region  with  varying  amount  of  noise  will  be 
interesting . 
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2.11  Boundary  Detection  in  Range  Pictures* 


S.  Inokuchi  and  R.  Nevatia 


Introduction 

Availability  of  range  i.e.  the  distances  of  observed  points  of  a 
scene  from  the  viewer,  simplifies  many  scene  analysis  problems 
including  segmentation  and  shape  analysis.  However,  extraction  of 
boundaries  of  object  surfaces  is  still  a  non-trivial  problem  and  many 
techniques  have  been  described  in  previous  work  [1-4].  Most  of  the 
techniques  may  be  said  to  use  a  "region  approach";  they  attempt  to 
isolate  planar  regions  or  grow  smooth,  curved  regions  from  smaller 
planar  regions.  Use  of  edge  like  processing  has  been  limited  to  the 
use  of  "jump  boundaries",  corresponding  to  a  large  range  discontinuity 
caused  by  large  separation  between  objects  [1,4]. 

In  this  paper,  we  describe  a  more  comprehensive  approach  to 
boundary  detection  using  edge  processing.  In  a  previous  comparison  of 
segmentation  of  intensity  images,  Nevatia  and  Price  conclude  that  the 
edge  and  region  techniques  are  complementary  and  each  is  suited  to 
particular  tasks  [5],  It  is  hoped  that  similar  observations  will  hold 
for  range  data  processing.  In  particular,  the  edge  methods  should 
work  better  with  thin  and  long  surfaces  thus  preserving  fine  detail, 
and  be  more  adaptable  to  curved  surfaces. 


*This  work  is  largely  supported  by  a  fellowship  to  Dr.  S.  Tnokuchi  , 
who  is  visiting  USC  from  Osaka  University,  Osaka,  Japan.  It  is 
included  here  due  to  possible  applications  in  the  DARPA  IU  work. 
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Boundary  Detection 

Our  boundary  detection  process  proceeds  in  the  following  steps. 

a .  Preliminary  Edge  and  Line  Detection 

The  first  stage  of  processing  uses  an  edge  and  line  detection 
technique  developed  for  intensity  images  and  described  in  T6].  The 
image  is  first  convolved  with  ideal  edge  masks  in  six  directions,  30 
degrees  apart,  the  maximum  output  determining  the  edge  magnitude  and 
the  direction.  The  edge  magnitude  and  directions  are  used  to  obtain 
thinned  and  thresholded  binary  edges.  In  our  experiments,  we  have 
used  edge  masks  suitable  for  detecting  step  edges.  Thus,  edges  along 
intersection  of  two  surfaces  may  not  be  detected  reliably.  Also, 
false  edges  m  ly  be  detected  on  surfaces  with  a  high  slope  from  the 
viewer.  The  latter  difficulty  could  be  avoided  by  use  of  second 
derivative  or  "bar"  masks.  We  expect  to  use  such  masks  for  future 
work.  The  remainder  of  the  technique  is  not  expected  to  be  affected 
by  this  change. 

The  binary  edges  are  linked  to  their  neighbors  with  similar  edge 
directions,  and  finally  the  linked  segments  are  approximated  by 
piecewise  linear  segments.  Figure  1(a)  shows  a  range  image  and  Fig. 
1(b)  the  corresponding  intensity  image.  Figure  2(a)  shows  the  edges 
detected  in  the  range  image  and  Fig.  2(b)  shows  the  resulting 
segments  (the  numeric  labels  in  Fig.  2(b)  are  explained  later). 

b.  Extension  of  Segments 

Fig.  2(b)  shows  many  segments  that  are  incomplete,  in  that  they 
terminate  without  connecting  to  other  segments.  In  usual  intensity 
image  processing,  proper  extension  of  such  segments  is  a  difficult 
task.  In  range  images,  we  can  use  3-D  information  to  facilitate  such 
extension.  The  extension  proceeds  in  three  steps,  corresponding  to 
the  three  types  of  "loose  ends"  labeled  in  Fig.  2(b). 
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i)  Loose  ends  marked  "1"  occur  near  an  intersection  of  two 
segments.  At  such  points,  the  local  edge  directions  are 

unreliable.  We  extend  such  segments  by  ignoring  edge  direction 
in  the  thinning  and  linking  steps  in  (a)  above,  and  requiring  a 
larger  threshold  on  the  edge  magnitudes  instead. 

ii)  Loose  ends  labeled  "2"  occur  near  "T“  junctions.  Again,  the 
connecting  of  two  segments  is  inhibited  in  step  (a)  due  to 
restrictions  on  local  edge  directions,  which  are  now  relaxed  for 
extension . 

iii)  Loose  ends  labeled  "3“  are  more  interesting  and  caused  by  the 
3-D  nature  of  range  data.  In  the  examples  shown,  the  detected 
segments  are  edges  between  a  surface  and  another  surface  far 
behind  it.  However,  along  the  extension  of  the  edge  segment, 
the  two  surfaces  come  nearer,  eventually  interesting  at  a 
corner,  ar.d  the  desired  extensions  are  not  detected  as  jump 
boundaries.  To  extend  such  segments,  we  compute  3-D  equations 
of  two  lines  parallel  to  the  segment,  but  on  either  sides  of  it 
(i.e.  lying  in  two  separate  planes).  If  these  two  lines 
intersect,  and  the  intersection  is  along  the  extrapolated 
boundary  segment,  and  there  exist  edge  pixels  whose  3-D  position 
is  near  the  computed  intersection,  then  the  boundary  segment  is 
extended . 

Figure  3  shows  the  results  of  all  the  above  three  extensions 
applied  to  segments  of  Fig.  2(b). 

c.  Radial  Line  Detection 

Some  of  the  interior  boundaries,  corresponding  to  intersection  of 
two  object  surfaces  may  be  still  missing.  Typically,  this  happens  for 
short  segments  at  the  corners.  Now,  for  each  corner,  we  examine  if  a 
line  can  be  found  in  some  direction  such  that  the  surfaces  on  its  two 


sides  correspond  to  different  planes.  A  window  (15x15)  is  centered 
around  each  corner  and  we  find  the  part  of  the  window  containing  the 
object  surfaces  (by  using  the  orientation  of  the  segments).  This 
region  is  divided  in  two  subregions  R^  and  R2  by  a  line  through  the 
corner  at  chosen  angles  (see  Fig.  4(a)).  For  each  choice  of  9  , 
planes  are  fit,  in  the  least  mean  squared  error  sense,  using  the  range 
data.  If  an  additional  edge  exists  at  this  corner,  for  some  angle,  9  , 
the  total  error  of  fitting  planes  to  R^  and  R^  should  be  small.  If 
the  r  itio  of  the  minimum  to  the  maximum  error,  for  different  9  ,  is 
less  than  1/2,  we  construct  a  line  segment  at  the  angle  corresponding 
to  the  minimum  error.  The  case  of  a  window  containing  two  corners  is 
shown  in  Fig.  4(b).  The  radial  line  detector  is  also  applied  to  loose 
ends  (an  example  is  the  segment  near  the  wheel  of  the  cart) .  Tn  this 
case,  the  detected  line  may  be  anywhere  in  the  window.  If  the  missing 
segment  is  long,  we  need  to  "follow"  the  newly  detected  segment. 
However,  this  last  step  is  not  necessary  for  our  example  and  we  have 
not  yet  implemented  such  a  line  following  procedure. 

Figure  5  shows  all  segments  resulting  after  application  of  the 
radial  line  detector.  Note  that  the  segments  detected  at  this  step 
are  undirected  (the  previous  segments  are  directed  in  that  the  object 
surfaces,  assumed  to  be  nearer  to  the  viewer  than  the  background,  are 
to  the  left  of  them) . 

Our  boundary  detection  technique  is  similar  in  spirit  to  the 
heterarchical  technique  of  Shirai  for  intensity  images  of  polyhedral 
objects  [7 ]  . 

Region  Descriptions 

Finding  boundaries  of  closed  regions  is  relatively  simple,  if 
complete  boundary  segments  are  found.  We  are  currently  finishing  the 
region  tracing  programs.  Interestingly,  some  region  properties  can  be 
inferred  by  observing  the  orientation  of  segments  along  a  region  as 
shown  in  Fig.  6.  Basically,  the  parts  and  spaces  can  be 
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Fig.  3.  Edges  after  extension  of 
line  segments. 


Fig.  4.  Radial  line  detection 
a)  and  b) . 


Fig.  5.  Edges  after  radial  line  detection. 
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differentiated  by  whether  the  region  is  traced  in  clockwise  or 
countr-clockwi se  direction.  Occlusion  is  indicated  by  not  all 
segments  pointing  in  the  same  direction,  compared  to  the  direction  of 
region  tracing.  Related  line  labeling  analysis  may  be  found  in  [81. 
Figure  7  shows  some  of  the  regions  in  our  example  that  can  be  inferred 
from  them.  This  part  of  our  process  is  currently  being  implemented. 
Eventually,  we  hope  to  better  describe  occluded  sub-parts,  and  be  able 
to  merge  such  sub-parts  when  appropriate. 

Conclusions 


We  have  described  a  technique  to  extract  fairly  complete 
boundaries  from  range  data,  exploiting  the  properties  of  such  data. 
Detailed  comparative  evaluations  with  region  based  techniques  are  yet 
to  be  made.  However,  we  expect  our  technique  to  be  complementary  and 
more  suitable  for  certain  tasks. 

References 


1.  D.  Nitzan,  A.E.  Brain  and  R.O  Duda,  "The  Measurement  and  Use  of 
Registered  Reflectance  and  Range  Data  in  Scene  Analysis,"  Proc.  IEEE , 
Vol .  69,  No.  2,  Feb.  77,  pp.  206-220. 

2.  R.O.  Duda,  D.  Nitzan  and  P.  Barrett,  "Use  of  Range  and  Reflectance 
Data  to  Find  Planar  Surface  Regions,"  IEEE  Trans ,  -PAMI ,  Vol.  1, 
No.  3,  July  1979,  pp.  259-271. 

3.  M.  Oshima  and  Y.  Shirai,  "A  Scene  Description  Method  Using 
Three-Dimensional  Information,"  Pattern  Recognition  Journal,  1979, 
pp.  9-17. 

4.  R.  Nevatia  and  T.O.  Binford,  "Description  and  Recognition  of 
Curved  Objects,"  Artificial  Intell igence,  Vol.  8,  No.  1,  Feb.  77, 
pp.  77-98. 


157 


5.  R.  Nevatia  and  K.  Price,  "Locating  Structures  in  Aerial  Images," 
Proc.  of  the  IJCPR,  Kyoto,  Japan,  Nov.  78  pp.  686-690. 

6.  R.  Nevatia  and  K.R.  Babu,  "Linear  Feature  Extraction  and 
Description,"  Proc.  of  IJCAI-79,  Tokyo,  Japan,  August  1979, 
pp.  639-641. 

7.  Y.  Shirai,  "Analyzing  Intensity  Arrays  Using  Knowledge  About 
Scenes,"  in  The  Psychology  of  Computer  Vision,  P.H.  Winston,  Editor, 
McGraw-Hill,  1975. 

8.  K.  Sugihara,  "Range  Data  Analysis  Guided  by  a  Junction 
Dictionary,"  Artificial  Intelligence  Journal,  Vol.  12,  1979, 
pp.  41-69. 


158 


3.  HARDWARE  IMPLEMENTATION  OF  IU  ALGORITHMS 


3.1  ADVANCED  IMAGE  UNDERSTANDING  USING  LSI  AND  VLSI 

S.D.  Fouse,  G.R.  Nudd,  and  V.S.  Wong 

Hughes  Research  Laboratories 
Malibu,  California  90265 


Abstract 

We  describe  here  the  work  undertaken  at  the  Hughes  Research 
Laboratories,  Malibu,  in  support  of  the  DARPA  Image  Understanding  (IU) 
program.  This  report  covers  the  period  from  October  1979  through 
April  1980,  and  as  such  represents  a  transition  period  for  us.  Our 
work  prior  to  this  period  has  been  concerned  with  the  investigation  and 
demonstration  of  large  scale  integration  (LSI)  microelectronic  technology 
for  image  understanding.  The  principal  aims  of  this  work  were  the 
design,  fabrication  of  and  demonstration  of  high-speed  primitives  for 
real-time  processing  with  low-level  operators.  This  work  is  continuing 
in  that  we  are  completing  a  performance  evaluation  of  the  fourteen  low- 
level  operators  already  implemented  on  the  program,  and  interacting  with 
other  groups  in  the  program,  (including  USC,  MIT,  and  Stanford)  is 
determined  how  these  might  perform  in  a  full-scale  IU  system.  The  second 
issue  which  is  now  becoming  the  mainstream  of  the  program  is  a  detached 
investigation  of  the  applicability  of  very  large  scale  integration  (VLSI, 
>50,000  gates/chip)  to  higher  level  systems.  The  aims  and  our  progress 
on  this  is  described  below. 
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I.  INTRODUCTION 

Our  work  this  period  on  the  image  understanding  (IU)  program  has 
been  divided  into  two  distinct  topics:  the  performance  investigation  of 
the  LSI  circuit  developed  for  the  low-lead  operators  and  an  investigation 
and  analysis  of  the  potential  of  VLSI  for  IU.  We  are  starting  a  detailed 
performance  review  to  characterize  the  throughput  accuracy  and  dynamic 
range  of  each  of  the  14  nrimitives  developed  to  date  on  the  program. 

This  will  enable  us  to  determine  how  they  might  interface  with  either  a 
host  machine  or  operate  as  part  of  a  more  complex  VLSI  system.  (To  this 
end,  we  are  supporting  both  USC  and  the  Artificial  Intelligence  Labora¬ 
tory  of  MIT  in  their  vision  projects.) 

The  second  topic  dealing  with  the  application  of  VLSI  to  IU  has 
involved  three  issues,  on  each  of  which  we  have  made  some  progress: 

(1)  a  directed  graph  analysis  of  candidate  systems,  (2)  an  investigation 
of  the  applicability  of  residue  operations  for  VLSI,  (3)  development  of 
a  simulation  capability  for  the  processors.  Initially,  we  have  chosen 
to  look  at  three  systems:  a  line  finder,  texture  analyzer,  and  a  seg¬ 
mentation  scheme.  Details  of  each  topic  are  given  below. 

II.  TEST  CHIP  III  PERFORMANCE  EVALUATION 

One  aspect  of  our  work  this  year  has  been  to  continue  the  testing 

and  evaluation  of  the  third  test  chip,  a  CCD-MOS  design.  This  chip,  as 

previously  described  in  Ref.  1,  has  five  functions:  a  5x5  programmable 

convolution,  a  7x7  mask  programmable  convolution  currently  implemented 

as  an  edge  detector,  a  3x3  Laplacian  operator,  a  5-element  sort  for 

median  filtering,  and  a  26x26  bipolar  convolution.  Preliminary  test 
2 

results  demonstrated  that  the  circuits  were  functionally  correct  (e.g., 
the  5x5  kernel  was  programmable,  the  sort  circuit  could  perform  the  sort, 
and  the  26x26  convolution  had  the  proper  impulse  response) .  Results 
included  images  that  were  processed  by  the  chip  at  pixel  rates  of 
^'20  kHz.  Our  current  work  has  been  to  extend  the  testing,  both  in  terms 
of  quantitative  performance  evaluation  and  in  terms  of  achievable  pixel 
rate. 
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Some  of  the  plans  we  have  for  this  comprehensive  testing  include: 

(1)  Calculation  of  dynamic  range  of  programmable  weights 

(2)  Determination  of  linearity  of  weights 

(3)  Experimental  analysis  of  sample-and-hold  circuit 

(4)  Dynamic  range  and  linearity  of  sort  circuit 

(5)  Performance  at  higher  speeds. 

To  help  complete  this  evaluation  in  a  timely  manner,  we  have  trained  a 
technician  to  operate  these  circuits.  This  should  enable  a  much  more 
thorough  job  to  be  done.  In  addition,  a  new  lot  of  chips  is  currently 
being  processed,  giving  us  a  larger  sample  of  the  circuits  to  work  with 

r 

as  well  as  possibly  improved  chips. 

III.  DEVELOPMENT  OF  INTERMEDIATE-LEVEL  VLSI  IMAGE-UNDERSTANDING 

SYSTEM 

A  major  goal  of  the  current  phase  of  our  work  on  hardware  for  IU 
systems  is  to  determine  the  impact  that  VLSI  technology  can  have  in  this 
area.  Our  method  for  accomplishing  this  is  shown  in  Figure  1.  We  have 
selected  candidate  IU  systems  that  extend  beyond  the  low-level  processing 
which  has  been  implemented  previously.  The  selected  systems  have  been 
characterized  using  directed  graphs,  which  help  to  make  explicit  the 
parallelism  that  is  inherent  in  the  algorithms.  The  graphs  will  also 
aid  in  determining  the  commonality  between  the  system,  thus  enabling  us 
to  select  a  partioning  in  which  functions  constructed  can  be  used,  as  is, 
in  several  systems. 

In  addition  to  our  directed-graph  analysis,  we  have  also  been 
working  in  two  areas  to  support  the  development  of  VLSI  architectures: 
residue  computation  techniques  and  the  acquisition  and  familiarization  of 
a  computer  simulation  program  (ECSS).  The  work  on  residue  computation 
as  well  as  future  work  in  defining  an  operating  environment  (e.g.,  will 
this  processor  be  used  as  a  peripheral  to  a  general-purpose  computer  or 
as  a  stand-alone  system),  will  he  fed  into  the  selection  of  candidate 
architectures.  Following  the  selectio..  j.  the  initial  architectures,  an 
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Figure  1.  Program  plan  for  development  of  intermediate-level  VLSI  image  understanding 
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iterative  process  of  simulation  using  ECSS,  performance  analysis,  and 
re-specification  of  the  architecture  will  be  undertaken.  Once  this 
process  converges  to  an  optimal  design,  we  will  take  the  further  step  of 
generating  a  fairly  detailed  design  of  a  chip  set  for  one  of  the 
algorithms. 

A  more  detailed  description  of  the  work  that  has  been  done  to  date 
on  this  plan  follows,  including  as  well  a  preliminary  sizing  of  one  of 
the  candidate  systems.  This  was  included  to  provide  a  little  more 
insight  into  the  scale  of  the  problem  we  are  trying  to  tackle. 

A.  Directed  Graph  Analysis  of  Candidate  Systems 

Three  systems  involving  intermediate  load  processing  were  selected 
for  the  purpose  of  developing  architectures  and  utilizing  commonality 
between  them.  The  systems  chosen  were  a  line-finding  system,  a  texture- 
classification  system,  and  a  segmentation  system.  To  be  able  to  perform 
the  directed  graph  analysis,  algorithms  must  be  specified;  the  following 
algorithms  were  chosen: 

(1)  Nevatia  line  finder 

(2)  Laws  texture  classification  system 

(3)  Ohlander /Price  segmenter. 

The  graph  analysis  was  carried  out  using  papers  describing  the  algorithms 

3  A  5 

and  by  talking  to  the  individuals  involved  in  their  development.  *  ’ 

The  graphs  for  the  three  systems  are  shown  in  Figures  2,  3,  and  4.  Con¬ 
structing  the  graphs  involves  identifying  logical  function  blocks  and 
the  data  inputs  and  outputs  of  each  block.  To  indicate  the  throughput 
required,  the  word  length  and  the  data  rates  are  indicated  on  the  arcs 
connecting  the  function  blocks.  For  example,  for  the  Nevatia  line  finder 
in  Figure  3,  the  input  data  is  an  eight-bit  parallel  word  and  the  data 
rate  is  NxN  pixels  per  frame  multiplied  by  F  frames  per  second. 

Once  the  graphs  are  constructed,  we  can  begin  to  see  the  parallelism 
potentials  of  the  algorithms  as  well  as  the  common  functions  between  the 
systems.  For  example,  for  the  Laws  texture  system,  it  is  apparent  that 
the  processing  from  the  5x5  convolutions  through  to  the  energy  measure 
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involves  M  independent  signal  paths,  indicating  that  a  logicalpartition 
would  be  to  put  the  5x5  convolution,  the  normalization,  and  the  large 
window  energy  measure  on  a  single  chip  (if  possible) ;  to  realize  this 
system,  we  need  only  make  M  copies  of  the  chip.  The  advantage  of  this 
type  of  partitioning  is  that  the  communication  to  and  from  the  chip  is 
minimized.  The  several  similarities  between  the  systems  are  made  obvi¬ 
ous  by  the  graphs.  There  are  parallel  convolutions  in  both  the  line- 
finder  and  the  texture  systems.  A  less  obvious  similarity  is  between 
the  edge-linker  function  in  the  line  finder  and  the  connected  region 
finder  in  the  segmentation  system.  Both  functions  play  the  same  role  of 
linking  data  points  into  a  single  element,  either  a  line  or  a  region. 

B.  Residue  Techniques  for  an  IU  Processor 

To  appreciate  our  interest  in  residue  arithmetic  and  its  applica¬ 
tion  to  IU, it  is  necessary  to  understand  the  mathematical  foundations. 

To  explain  it,  some  mathematical  notation  must  be  introduced.  The 
expression  "A  Mod  M"  is  equal  to  the  remainder  when  A  is  divided  by  M. 
Both  A  and  M  must  be  integers  and  A  mod  M  must  also  be  an  integer  between 
0  and  M-l.  Using  modern  algebra,  it  is  possible  to  show  two  facts. 

First  is  the  expression 

((A  mod  M)  +  (B  mod  M) )  mod  M  =  (A  +  B)  mod  M  . 

This  implies  that  we  can  calculate  the  remainder  of  a  sum  by  adding  the 
remainders  themselves.  The  second  fact  is  that  the  k-tuple  (A  mod  Ml, 

A  mod  M2,..., A  mod  Mk)  is  a  unique  representation  of  A  if  two  conditions 
are  met:  the  M^'s  are  relatively  prime,  and  A  is  less  than  the  product 
of  the  M^'s.  These  two  conditions  allow  two  large  numbers  to  be  added 
by  adding  several  pairs  of  smaller  numbers,  the  smaller  numbers  being 
remainders.  This  can  also  be  said  for  multiplication. 

To  perform  residue  operations,  three  steps  are  involved.  First, 
the  binary  operands  must  be  converted  to  remainders  of  several  bases 
(M^'s).  Next,  these  remainders  must  be  added  in  a  modular  fashion 
(i.e.,  the  sum  of  the  two  remainders  must  also  be  a  remainder,  which  is 
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to  say  that  one  must  take  the  remainder  of  the  sum).  Finally,  the  result 
must  be  converted  back  to  binary,  which  seems  to  be  the  most  difficult 
of  the  three  steps.  In  addition  to  this  overhead  of  conversion,  it  is 
difficult  to  make  decisions  while  in  the  residue  representation  since, 
unlike  binary,  there  is  no  ordering  associated  with  the  representation. 

The  overhead  costs  associated  with  the  residue  technique  for  com¬ 
putation  makes  it  inappropriate  for  many  applications.  But  often  it  is 
easier  to  do  several  operations  with  small  operands  than  to  do  one  oper¬ 
ation  with  large  operands.  So  if  enough  operations  can  be  done  between 
the  two  conversion  steps,  the  process  as  a  whole  could  be  cheaper  to  do 
than  the  equivalent  process  using  binary  arithmetic.  For  IU  work,  the 
bulk  of  the  computation  is  done  at  the  low  level  and  consists  mostly  of 
signal  processing  requiring  no  decisions  to  be  made.  For  algorithms 
such  as  a  line  finder  that  require  several  parallel  convolutions,  there 
are  numerous  operations  that  can  be  performed  in  the  residue  form. 
Additionally,  we  are  dealing  with  integer  quantities,  which  are  ideally 
suited  to  residue  techniques.  In  terms  of  building  hardware,  a  residue 
processor  would  have  relatively  simple  logic  performing  the  bulk  of  the 
operations  and  the  complex  logic  would  be  concentrated  at  the  input  and 
the  output  of  the  system.  This  is  an  advantage  in  terms  of  reliability 
and  testability. 

Our  goal  for  this  program  is  to  examine  the  different  systems  that 
we  have  chosen  and  to  determine  the  feasibility  of  using  the  residue 
technique.  The  study  will  be  done  in  a  parametric  fashion,  so  that,  by 
using  the  parameters  of  a  particular  system,  one  can  determine  if  residue 
will  be  an  advantage.  Areas  that  we  have  already  begun  to  look  at  or 
will  be  looking  at  include  conversion  logic,  dynamic  range  requirements 
(rounding  in  residue  makes  little  sense),  logic  for  adding  and  multiply¬ 
ing  in  residue,  methods  for  scaling  intermediate  results,  and  methods  of 
representing  the  numbers  in  residue  to  facilitate  the  logic  for  calculations. 

C.  ECSS  Simulation 

As  the  density  of  micro-electronic  circuits  increases,  the  systems 
that  can  be  built  on  these  chips  may  become  very  sophisticated.  To 


optimize  the  design  and  to  be  able  to  predict  performance,  it  is  desirable 
to  simulate  the  operation  before  the  circuits  are  actually  constructed  in 
hardware.  For  our  work  on  systems  design  at  HRL,  we  have  chosen  to  use 
a  program  developed  specifically  for  computer  system  simulation  at  the 
RAND  Corporation  called  ECSS  (Extendable  Computer  Simulation  System). 

ECSS  is  written  in  a  natural,  English-like  language,  and  it  is  both 
descriptive  and  flexible  to  use.  It  stresses  a  procedural  format,  and 
the  syntax  of  the  language  is  constructed  so  that  it  can  compactly  define 
and  specify  the  components  of  a  computer  system  (storage  devices,  I/O 
devices,  CPUs,  printers)  as  well  as  interactions  among  the  components 
(job  scheduling,  resource  allocation,  interprocess  communication).  To 
indicate  performance  and  bottlenecks  of  the  computer  system  being  modeled, 
ECSS  outputs  statistics  showing  the  jobs  activated,  percent  utilization 
of  the  devices,  and  length  of  queues  for  the  various  resources. 

The  main  advantages  of  using  ECSS  over  other  simulation  programs 
include: 

•  Flexibility.  ECSS  provides  a  variety  of  techniques  for 
modeling  systems.  The  user,  depending  on  his  interests, 

may  focus  on  modeling  program  behavior,  system  I/O  transfers, 
or  arbitrary  activities  running  on  a  particular  device. 

•  Extendability.  ECSS  is  built  on  SIMSCRIPT  II,  a  program 
language  designed  for  discrete  event  simulation.  There¬ 
fore,  ECSS  has  all  the  provisions  of  data  types  and 
procedural  statements  of  SIMSCRIPT.  This  allows  the  user 
to  extend  definitions  of  ECSS  structures  (devices,  jobs, 
transmissions)  and  add  new  commands  as  desired. 

•  Programmability.  When  the  models  provided  by  ECSS  do 
not  meet  the  requirements  of  the  user,  it  is  possible 

to  modify  the  source  code  for  ECSS.  The  service  routines 
are  written  in  modular  fashion  in  SIMSCRIPT,  thus 
facilitating  modification. 

At  HRL,  we  are  presently  using  ECSS  to  simulate  computer  architec¬ 
tures  for  our  IU  program  and  to  investigate  specialized  VLSI  designs  for 
a  line  finder,  a  texture  system,  and  a  segmentation  system.  We  antici¬ 
pate  that  ECSS  will  be  a  very  useful  tool  for  developing  the  architec¬ 
tures,  predicting  the  performance  of  the  systems,  and  in  general  for 
communicating  and  documenting  the  designs.  We  also  have  interest  in  the 


area  of  specialized  cellular  array  processors  for  executing  concurrent 
algorithms.  In  this  respect,  we  believe  that  ECSS  would  be  very  useful 
in  the  design  of  such  an  array  processor. 

D.  Preliminary  Chip  Sizing  of  Line-Finder  System 

To  obtain  a  clearer  idea  of  the  magnitude  of  the  problem  we  are 
trying  to  solve,  as  well  as  to  determine  the  pacing  items  for  fabricating 
the  system,  we  have  decided  to  determine  roughly  how  many  gates  the  line- 
finder  system  would  need.  The  gate  count  was  divided  into  gates  for 
memory  and  gates  for  random  logic.  We  then  used  these  figures,  along 
with  a  state-of-the-art  number  for  gates  per  chip  for  both  random  logic 
and  memory,  to  determine  an  approximate  total  number  of  chips  for  the 
system.  Additionally,  we  calculated  a  chip  count  for  the  predicted 
1985-1990  technology. 

The  gate  count  for  the  line  finder  was  done  for  each  block  on  the 
directed  graph.  These  functions  include: 


(1) 

Kernel  generation 

(2) 

5x5  convolutions 

(3) 

Edge  detection 

(4) 

Edge  thinning 

(5) 

Predesessor  and  successor  generation 

(6) 

P  and  S  memory 

(7) 

Edge  linking. 

Rough  designs  were  performed  for  the  convolutions,  edge  detection,  and 

edge  thinning  to  estimate  the  number  of  gates.  The  P  and  S  generation 

can  be  estimated  by  replicating  the  edge-thinning  logic  six  to  ten  times 

since  it  is  making  similar  comparisons  as  for  edge  thinning  but  more  of 

them.  The  logic  for  edge  linking  gets  very  involved,  since  it  must  be 
able  to  generate  addresses  for  random  access  of  the  P  and  S  memory.  In 
addition,  logic  for  each  pass  must  be  generated,  since  each  pass  performs 
a  different  function.  So  for  the  purposes  of  this  sizing,  the  edge 
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linking  was  not  estimated.  Table  1  shows  the  gate  counts  for  each  of 
the  functions.  Notice  that,  for  the  convolutions,  we  sized  two  types, 
one  with  no  multiplies  and  one  with  six  multiplies.  This  is  because  the 
algorithm  we  are  sizing  has  six  convolutions,  two  which  can  be  scaled 
such  that  all  the  weights  are  either  zero  or  one,  and  four  which  can  be 
scaled  such  that  all  but  six  of  the  weights  are  either  zero  or  one. 

The  current  state  of  the  art  for  MOS  technology  allows  us  to  put 
64K  bits  of  memory  or  20,000  gates  of  random  logic  on  a  single  chip. 

Using  these  figures,  we  are  able  to  partition  the  algorithm  onto  a  set 
of  chips  and  thus  calculate  a  total  number  of  chips  for  the  system.  One 
possible  partitioning  is  as  follows: 

CHIP  FUNCTIONS 

1  5  line  kernel,  2  convolutions  (no  multiply) 

2-5  One  5x5  convolution  (6  multiplies) 

6  Edge  detection,  3  line  kernel,  thinning  logic 

7  3  line  kernel,  P  and  S  generation 

8-27  P  and  S  memory  (40  64K  chips) 

48  Edge  linking  logic. 

It  is  fairly  obvious  to  see  that  the  two  pacing  items  (not  counting  the 
edge-linking  logic,  which  we  have  not  looked  at)  are  the  P  and  S  memory 
and  the  convolution  calculations.  The  convolution  complexity  can  be 
simplified  by  reducing  the  accuracy  required  on  the  weights  of  the  kernel. 
This  would  mean  that  a  logic  multiplier  could  be  fast  enough,  and  the 
memory  for  look  up  table  multiplies  would  not  be  needed.  The  require¬ 
ment  for  the  P  and  S  memory  is  embedded  in  the  algorithm,  and  thus  the 
algorithm  would  need  to  be  altered  if  we  wanted  to  reduce  this  memory. 

We  can  perform  this  same  analysis  assuming  some  future  technology. 

We  expect  that  in  several  years  we  will  be  able  to  achieve  1  Mbit  of 
memory  or  0.5  million  gates  of  random  logic  on  a  single  chip.  This 
assumes  1-um  feature  size  as  compared  to  5-um  features  for  current  tech¬ 
nology.  For  this  case,  we  get  the  following  partitioning: 
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CHIP  FUNCTION 

1  All  random  logic  and  kernel  generation  memory 

2-4  P  and  S  memory. 

Here  we  can  have  8  lines  for  the  input,  8  lines  for  the  output,  and  some 
address  and  data  lines  for  the  P  and  S  memory,  avoiding  any  pin-out 
problems. 

The  papers  presented  as  Appendices  A  and  B  were  prepared  during 
this  quarter. 


Table  1.  Gate  Count  for  Line  Finder 


Function 

Gate  Count 

Number  Used 

Total 

5  line  kernel 

4x5121x8  =  16K 
bits 

1 

16K  bits  memory 

5x5  convolution 
no  multiplies 

1.25K 

2 

2.5K 

5x5  convolution 

6  multiplies2 

13.25K 

4 

53K 

Edge  detection 

0.35K 

1 

0.35K 

3  line  kernel 

2x512x12  =  6K 
bits 

2 

12K  bits  memory 

Edge  thinning 

0.35K 

1 

0.35K 

P  and  S 
generation 

3K 

1 

3K 

P  and  S 
memory 

512x512x10 
=  2.5  Mbit 

1 

2.5  Mbit  memory 

1  -  Assumed  image  size  -  512x512 

2  -  Multiplies  accomplished  using 

8x256  bit  ROM. 

The  following  papers  describing  this  work  were  prepared  during  the  last 

■is  months. 


s.D.  rouse,  G.R.  Nudd  and  P.A.  Nygaard,  "Implementation  of  Image 
Preprocessing  Functions  Using  CCD  LSI  Circuits,"  Proc.  SPIE  Tech.  Symp.  -  IR 
Sensor  Technology,  Vol.  225,  Washington,  D.C.,  April  1980. 

-  ’.udd.  "Image  Understanding  Architectures,"  National  Computer  Conference, 
.  ■  iTi,  Ca .  ,  May  1980. 


172 


REFERENCES 


1.  G.R.  Ntidd,  P.A.  Nygaard,  S.D.  Fouse,  and  T.A.  Nussmeier, 
"Implementation  of  Advanced  Real-Time  Image  Understanding 
Algorithms,"  Proceedings  Image  Understanding  Workshop,  DARPA, 

April  1979. 

2.  G.R.  Nudd,  S.D.  Fouse,  T.A.  Nussmeier,  and  P.A.  Nygaard, 

"Development  of  Custom-Designed  Integrated  Circuits  for  Image 
Understanding,"  Proceedings  of  Image  Understanding  Workshop,  DARPA, 
November  1979. 

3.  R.  Nevatia  and  K.  Ramesh  Babu,  "An  Edge  Detection,  Linking,  and  Line 
Finding  Routine,"  USC  IPI  Semiannual  Report,  September  1978. 

4.  K.I.  Laws,  "Textured  Image  Segmentation,"  Ph.D.  Thesis,  USC, 

January  1980. 

5.  R.  Ohlander,  K.  Price,  and  D.  Raj  Reddy,  "Picture  Segmentation 
Using  a  Recursive  Region  Splitting  Method,"  Computer  Graphics  and 
Image  Processing,  1978. 


173 


4.  RECENT  INSTITUTE  PERSONNEL  PUBLICATIONS  AND  PRESENTATIONS 

1.  A.  Armand ,  D.  Boswell,  A. A.  Sawchuk,  B.H.  Soffer,  T.C.  Strand, 
"Real-Time  Parallel  Optical  Analog-to-Dig ital  Conversion,"  Opt i cs 
Letters ,  Vol.  5,  March  1980. 

2.  A.  Armand,  A. A.  Sawchuk,  T.C.  Strand,  "Nonlinear  Optical 

Processing  with  Halftones:  Accurate  Predictions  for  Degradation 
and  Compensation ,"  submitted  to  Appl ied  Opt i cs . 

3.  A.  Armand,  A. A.  Sawchuk,  T.C.  Strand,  "Real-Time  Paralle"1 

Logarithmic  Filtering,"  submitted  to  Appl i ed  Opt i cs . 

4.  J.  Bescos,  I.  Glaser,  A. A.  Sawchuk,  "Restoration  of  Color  Images 
Degraded  by  Chromatic  Aberrations,"  submitted  to  Appl ied  Opt i cs , 
Vol.  19. 

5.  D.  Boswell,  P.  Chavel  ,  A.M.  Lackner,  A. A.  Sawchuk,  B.H.  Soffer, 

T.C.  Strand,  A . R .  Tanguay,  Jr.,  "Optical  Computing  with 

Variable-Grating  Mode  Liquid  Crystal  Light  Valves,"  Proceedings 
1980  International  Opt i ca 1  Computing  Conference/SPIE  P roceed i ngs , 
Vol.  232,  Washington,  D.C. ,  April  1980. 

6.  D.  Boswell,  P.  Chavel,  A. A.  Sawchuk,  B.H.  Soffer,  T.C.  Strand, 
A.R.  Tanguay,  Jr.,  "Parallel  Optical  Anal og-to-Di g  ital  Conversion 
Using  a  Liquid  Crystal  Light  Valve,"  Workshop  on  High  Speed  A/D 
Conversion,  Portland,  Oregon,  February  1980. 

7.  D.  Boswell,  A.M.  Lackner,  A. A.  Sawchuk,  B.H.  Soffer,  T.C.  Strand, 
A.R.  Tanguay,  Jr..  "Variable  Grating  Mode  Liquid  Crystal  Device 
for  Optical  Processing  and  Computing,"  Eighth  International 
Liquid  Crystal  Conference,  Kyoto,  Japan,  June-July  1980. 


174 


8.  D.  Boswell,  A. A.  Sawchuk,  B.H.  Soffer,  T.C.  Strand,  A.R.  Tanguay, 
Jr.,  "Variable  Grating  Mode  Liquid  Crystal  Device  for  optica] 
Processing,"  Proceedings  Society  of  Photo-0ptica3  Instrumentation 
Engineers  Los  Angeles  Technical  Symposium  -  Devices  and  Systems 
for  Optical  Signal  Processing ,  Vol .  218,  Los  Angeles,  February 
1980. 

9.  O.D.  Faugeras,  "An  Overview  of  Probabilistic  Relaxation  Theory 
and  Applications,"  P roceedings  of  the  NATO  ARI ,  D.  Reide] 
Publishing  Company,  June  1980. 

10.  O.D.  Faugeras,  "Autoregressive  Modeling  with  Conditional 
Expectations  for  Texture  Synthesis,"  submitted  to  the  5th  ICPR, 
March  1980. 

11.  O.D.  Faugeras,  "Decomposition  and  Decentralization  Techniques  in 
Relaxation  Labeling,"  submitted  to  Computer  Graphics  and  Image 
Processing ,  March  1980. 

12.  O.D.  Faugeras,  "Improving  Consistency  and  Reducing  Ambiguity  in 
Stochastic  Labeling:  An  Optimization  Approach,"  to  be  published 
in  the  IEEE  PAMI  Transactions ,  1980  (with  M.  Berthod) . 

13.  O.D.  Faugeras,  "Scene  Labeling:  An  Optimization  Approach,"  to  be 
published  in  Pattern  Recognition,  1980  (with  M.  Berthod). 

14.  O.D.  Faugeras,  "Using  Context  in  the  Global  Recognition  of 
of  Objects:  An  Optimization  Approach,"  to  be  published 
Proceedings  of  the  8 th  World  Computer  Congress  (IFIP  80), 

M.  Berthod) . 

15.  O.D.  Faugeras,  D.D.  Garber,  "Algebraic  Reconstruction  Techniques 
for  Texture  Synthesis,"  submitted  to  the  5th  ICPR,  March  19P0. 


a  Ret 
in  the 
( wi  th 


175 


16.  O.D.  Faugeras,  W.K.  Pratt,  "Decorrelation  Methods  of  Texture 
Feature  Extraction,"  to  be  published  in  the  IEEE  PAMT 
Transactions,  3980. 


17.  O.D.  Faugeras,  K.E.  Price,  "Semantic  Description  of  Aerial  Images 
with  Stochastic  Labeling,"  submitted  to  the  5th  ICPR,  March  1980. 


18.  L.M.  Frantz,  A. A.  Sawchuk,  W.  Von 
Measurement  in  Real  Time," 
pp.  3301-3306,  October  1979. 


der  Ohe,  "Optical  Phase 
Appl ied  Optics ,  Vol .  18, 


19.  S.  Inokuchi  ,  and  R.  Nevatia,  "Boundary  Detection  in  Range  Data," 
submitted  for  conference  presentation. 

20.  C.M.  Lo  and  A. A.  Sawchuk,  "Nonlinear  Restoration  of  Blurred 
Images  with  Poisson  Noise,"  1979  Annual  Meeting,  Optical  Society 
of  America,  Vol.  ^9,  pp.  1441-1442,  Rochester,  New  York,  October 
1979. 

21.  C.M.  Lo ,  A. A.  Sawchuk,  "Restoration  with  Poisson  Noise," 
Proceedings  of  the  IEEE ,  Vol.  68,  1980  (invited  paper  to  appear). 

22.  J.  Mantock,  A. A.  Sawchuk,  T.C.  Strand,  "An  Optical  Processor 
Applied  to  Cloud  Classification,"  Proceedings  Soci ety  of 
Photo-Optical  Instrumentat ion  Engineers  Los  Angeles  Techn i ca 1 
Sympos i urn  -  Devices  and  Systems  for  Opti cal  S ignal  Process i ng , 
Vol.  218,  Los  Angeles,  February  1980. 

23.  J.  Mantock,  A. A.  Sawchuk,  T.C.  Strand,  "Hybrid  Optical -Dig i tal 
Texture  Analysis,"  First  ASSP  Workshop  on  Two-Dimensional  Digital 
Signal  Processing,  Berkeley,  California,  October  1979. 


24.  J.  Mantock,  A. A.  Sawchuk,  T.C.  Strand,  "Hybrid  Opt i cal /Di g i tal 
Texture  Analysis,"  Optical  Engineering,  Vol.  39,  March/April 
1980,  (invited  paper). 


176 


•mm 


RRHM 


25.  J.  Michaelson,  A. A.  Sawchuk,  "Nonlinear  Optical  Processing  Using 
Liquid  Crystal  Light  Valves,"  Proceedings  Society  of 
Photo-Optical  Instrumentation  Engineers  Los  Angel es  Technical 
Symposium  -  Devices  and  Systems  for  Optical  S igna 1  Processi ng , 
Vol.  218,  Los  Angeles,  February  1980. 

26.  R.  Nevatia,  "Image  Understanding  Research  at  USC,"  Seminar 
presented  at  University  of  California,  Irvin,  Ca . ,  March  1980. 

27.  R.  Nevatia,  "Matching  of  Natural  Terrain  Scenes,"  with  C.  Clark, 
D.  Conti,  W.  Eckhardt,  T.  McCullogh  and  D.  Tseng,  submitted  for 
conference  presentation. 

28.  R.  Nevatia  and  K.R.  Babu,  "Linear  Feature  Extraction  and 
Description,"  accepted  for  publication  in  Computer  G raph i cs  and 
Image  Processing. 

29.  R.  Nevatia,  and  K.E.  Price,  "Locating  Structures  in  Aerial 
Images,"  submitted  for  Journal  publication. 

30.  R.  Nevatia  and  K.E.  Price,  "Symbolic  Representation  in  USC  IU 
System,"  Proc.  of  ARPA  Image  Understanding  Workshop,  Los  Angeles, 
Ca . ,  November  1979. 

31.  R.  Nevatia  and  A. A.  Sawchuk,  "Image  Processing  and  Understand ing 
Research,"  Seminar  at  Battel le  Memorial  Institute,  Richland, 
Wash .  ,  Ma rch  1980. 

32.  R.  Nevati-'  and  A. A.  Sawchuk,  "Progress  in  Image  Understanding 
Research  at  USC,"  Proceed i ngs  Image  Understand i ng  Workshop, 
pp.  149-151,  Los  Angeles,  Ca . ,  November  1  979. 

33.  F.  Vilnrotter,  R.  Nevatia  and  K.E.  Price,  "Structural  Description 
of  Natural  Textures,"  submitted  for  conference  presentation. 


177 


5.  RECENT  PH.D.  DISSERTATIONS 


The  abstracts  of  recent  Ph.D.  Dissertations  completed  with  DARPA 
support  are  listed  in  this  section.  The  dissertations  are  published 
as  USCIPI  technical  reports. 


USCIPI  Report  940,  January  1980 

The  problem  of  image  texture  analysis  is  introduced,  and  existing 
approaches  are  surveyed.  An  empirical  evaluation  method  is  appl ied  to 
two  texture  measurement  systems,  co-occurrence  statistics  and 
augmented  correlation  statistics.  An  "spatial -statistical "  class  of 
texture  measures  is  then  defined  and  evaluated.  it  leads  to  a  simple 
class  of  "texture  energy"  transforms,  which  perform  better  than  any  of 
the  preceding  methods.  These  transforms  are  very  fast,  and  can  be 
made  invariant  to  changes  in  luminance,  contrast,  and  rotation  without 
histogram  equalization  or  other  preprocessing. 

Texture  energy  is  measured  by  filtering  with  small  masks, 
typically  5x5,  then  with  a  moving-window  average  of  the  absolute  image 
values.  This  method,  similar  to  human  visual  processing,  is 
appropriate  for  textures  with  short  coherent  length  or  correlation 
distance.  The  filter  masks  are  integer-val  ued  and  separate,  and  can 
be  implemented  with  one-dimensional  or  3x?  convolutions.  The 
averaging  operation  is  also  very  fast,  with  computing  time  independent 
of  window  size. 
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Texture  energy  planes  may  be  linearly  combined  to  form  a  smaller 
number  of  discriminant  planes.  These  principal  component  planes  seem 
to  represent  natural  texture  dimensions,  and  to  be  more  reliable 
textur e  measures  than  the  texture  energy  planes. 

Texture  segmentation  or  classification  may  be  accomplished  using 
either  texture  energy  or  principal  component  planes  as  input.  This 
study  classified  15x15  blocks  of  eight  natural  textures.  Accuracies 
of  72%  were  achieved  with  co-occurrence  statistics,  r<5%  with  augmented 
correlation  statistics,  and  94%  with  texture  energy  statistics. 


5.2  Design  of  SVD/SGK  Convolution  Filters  for  Image  Processing 
Sang  Uk  Lee 


USCIPI  Report  950,  January  1980 


This  dissertation  describes  a  special-purpose  signal  processor 
for  performing  two-dimensional  convolution  with  a  minimum  amount  of 
hardware  using  the  concepts  of  singular  value  decomposition  (SVD)  and 
small  generating  kernel  (SGK )  convolution.  The  SVT  of  an  impulse 
response  of  a  two-dimensional  finite  impulse  response  (FIR)  filter  is 
employed  to  decompose  a  filter  into  a  sum  of  two-dimensional  separable 
linear  operators.  These  linear  operators  are  themselves  decomposed 
into  a  sequence  of  small  kernel  convolution  operators.  The  FVD 
expansion  can  be  truncated  to  a  relatively  few  terms  without 
significantly  affecting  the  filter  output. 

A  statistical  analysis  of  finite  word-length  effects  in  SVD/RGK 
convolution  is  presented.  Two  important  issues,  related  to  the 
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implementation  of  the  filters  in  cascade  form,  scaling  and  section 
ordering,  are  also  considered. 

Computer  simulation  of  image  convolution  indicates  that  12  bits 
are  required  for  the  SVD/SGK  accumulator  memory  and  1A  bits  are 
required  for  quantization  of  filter  coefficients  to  obtain  results 
visually  indistinguishable  from  full  precision  computation.  A 
normalized  mean  square  error  between  the  SVD/SGK  processed  output  and 
the  direct  processed  output,  is  chosen  as  an  objective  criterion 
function.  It  is  shown  that  a  subjective  visual  improvement  is 
obtained  by  resetting  the  output  mean  to  be  equal  to  the  input  mean. 

The  transformation  technique  developed  for  the  one-dimens ional 
case  is  used  to  parametrically  modify  the  cutoff  frequency  of  a 
baseline  SVD/SGK  convolution  filter.  A  detailed  discussion  of  the 
one-dimensional  case  is  presented,  and  its  appl icabi 1 i ty  to  SVD/SGK 
convolution  filters  is  described. 
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